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Abstract: IceTop, the surface component of the IceCube Neutrino Observatory at the South Pole, is an air shower array 
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the energy range from about 10^* eV to 10^* eV by exploiting the correlation between the shower energy measured in 
IceTop and the energy deposited by muons in the deep ice. Construction of IceCube, including the IceTop component, 
was completed in December 2010. The final detector configuration, first operation and performance experiences, the 
development of an analysis framework, and first results will be reported. 
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1 Introduction 

The Neutrino Observatory IceCube is a 1-km'^ detector sit- 
uated in the ice of the geographic South Pole at a depth 
of about 2000 m. IceTop, the surface component of Ice- 
Cube, is an air shower array covering an area of Ikm^. 
The prime purpose of IceTop is the determination of the 
mass composition of primary cosmic rays in the energy 
range from about 10^'' eV to 10^^ eV. In the 'knee' region, 
at several PeV, the spectral index of the observed cosmic 
ray energy spectrum changes. Several experiments found 
this change to be accompanied by a change in the chemi- 
cal composition of the primaries. However, details of the 
features are not well known. In particular, there may be 
at the high end of the IceTop energy range another change 
of the spectral index and an accompanying change of the 
composition, possibly indicating the transition from galac- 
tic to extra-galactic origin of cosmic rays. An improve- 
ment of the experimental situation in this energy range - 
between direct measurements with balloons and satellites 
and the highest energies tackled by experiments Uke HiRes 
and Auger - is one of the main goals of cosmic ray physics 
with IceCube. 

The mass determination from extended air showers (EAS) 
is notoriously difficult because the measurements are indi- 
rect and have to rely on models for the hadronisation pro- 
cesses. Observables sensitive to the primary mass compo- 
sition are mainly the height of the shower maximum (mea- 
sured through fluorescence, Cherenkov or radio emission) 
and the number of muons in a shower Conceming the 
muon rate, the highest energy muons stemming from the 



first interactions in the higher atmosphere are most closely 
correlated to the mass of the primary nucleus. IceCube, in 
combination with IceTop, offers the unique possibility to 
observe these muons, typically with initial energies above 
about 500 GeV, in the deep ice in coincidence with the 
mostly electromagneticaUy deposited shower energy mea- 
sured at the surface. This provides an exceptionally power- 
ful method for the determination of the mass composition. 

To scrutinize the dependence on hadronisation models, sev- 
eral alternative methods for studying mass composition 
have been developed by the IceCube collaboration. Other 
mass sensitive observables are for example: the shower ab- 
sorption in the atmosphere at different zenith angles, the 
number of dominantly low-energy muons in the surface de- 
tector, and other shower properties such as shower age and 
shower front curvature. 

The IceTop array has additionally been used to study high- 
Pt muons, PeV-gammas and transient events, such as the 
radiation effects of solar flares. It also serves as a veto for 
the detection of downward-going neutrinos with IceCube 
and for direction calibration. 

2 The detector 

The IceCube construction was completed in December 
2010. The results presented here are based on data taken 
with smaller detector configurations. 

IceCube: The main coinponent of Icecube is an array of 
86 strings equipped with 5160 light detectors in a volume 
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Figure 1: The IceCube Observatory with its components 
DeepCore and IceTop. 

of 1 km-^ at a depth between 1450 m and 2450 m (Fig. 1). In 
the lower part of the detector a section called DeepCore is 
more densely instrumented. The main purpose of IceCube 
is the detection of high energy neutrinos from astrophysical 
sources via the Cherenkov light of charged particles gener- 
ated in neutrino interactions in the ice or the rock below the 



IceTop: The IceTop air shower array is located above 
IceCube at a height of 2832 m above sea level, correspond- 
ing to an atmospheric depth of about 680 g/cm^. It consists 
of 162 ice Cherenkov tanks, placed at 81 stations and dis- 
tributed over an area of 1 km^ on a grid with mean spacing 
of 125 m (Fig. 1). In the center of the array, three stations 
have been installed at intermediate positions. Together 
with the neighbouring stations they form an in-fill array for 
denser shower sampling. Each station comprises two cylin- 
drical tanks, 10 m apart from each other, with a diameter of 
1.86 m and filled with 90 cm ice. The tanks are embed- 
ded into the snow so that their top surface is level with the 
surrounding snow to minimize temperature variations and 
snow accumulation caused by wind drift. However, snow 
accumulation (mainly due to irregular snow surfaces) can- 
not be completely avoided so that the snow height has to 
be monitored (see ref. [1]) and taken into account in simu- 
lation and reconstruction (currently this is still a source of 
non-negligible systematic uncertainties). 

Each tank is equipped with two 'Digital Optical Mod- 
ules' (DOMs), each containing a 10" photo multiplier tube 
(PMT) to record the Cherenkov light of charged particles 
that penetrate the tank. In addition, a DOM houses complex 
electronic circuitry supplying signal digitisation, readout, 
triggering, calibration, data transfer and various control 
functions. The most important feature of the DOM elec- 
tronics is the recording of the analog waveforms in 3.3 ns 
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Figure 2: Reconstruction of shower parameters from the 
lateral distribution. 



wide bins for a duration of 422 ns. DOMs, electronics and 
readout scheme are the same as for the in-ice detector. 

The two DOMs in each tank are operated at different PMT 
gains (1 • 10^ and 5 • 10^) to cover a dynamic range of more 
than 10^. The measured charges are expressed in units of 
'vertical equivalent muons' (VEM) determined by calibrat- 
ing each DOM with muons (see ref. [1]). 

To initiate the readout of DOMs, a local coincidence of 
the two high gain DOMs of a station is required. This re- 
sults in a station trigger rate of about 30 Hz compared to 
about 1600 Hz of a single high gain DOM at a threshold 
of about 0.1 VEM. The data are written to a permanent 
storage medium, and are thus available for analysis, if the 
readouts of six or more DOMs are launched by a local coin- 
cidence. This leads to a trigger threshold of about 300 TeV. 
Additionally, IceTop is always read out in case of a trigger 
issued by another detector component (and vice versa). For 
each single tank above threshold, even without a local co- 
incidence, condensed data containing integrated charge and 
time stamp are transmitted. These so-called SLC hits (SLC 
= 'soft local coincidence') are useful for detecting single 
muons in showers where the electromagnetic component 
has been absorbed (low energies, outer region of showers, 
inclined showers). 

For monitoring transient events via rate variations, the time 
of single hits in different tanks with various thresholds are 
histogrammed. 

3 Shower reconstruction 

For each triggered tank in an event, time and charge of 
the signal are evaluated for further processing. Likelihood 
maximisation methods are used to reconstruct location, di- 
rection and size of the recorded showers. In general, signal 
times contain the direction information, and the charge dis- 
tribution is connected to shower size and core location. The 
standard analysis requires five or more triggered stations 
leading to a reconstruction threshold of about 500 TeV. A 
constant efficiency is reached at about 1 PeV, depending 
on shower inclination. For small showers an effort was 
launched to decrease the threshold to about 100 TeV with 
a modified reconstruction requiring only three stations. 
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The lateral signal distribution is fitted by a function which 
describes the logarithm of the tank signals as a second or- 
der polynomial in the logarithm of the distance from the 
shower axis (Fig. 2). Characterizing the shower size the 
signal 5*125 at a reference radius R^ci — 125 m, coincid- 
ing with the grid spacing, has been used for the analyses 
presented in this paper Studies of alternative lateral distri- 
bution functions are reported in a separate contribution to 
this conference [2] . 

The true energy spectrum is obtained by unfolding the 5*125 
distribution in different zenith angular ranges. Since the 
unfolding matrices depend on the primary mass composi- 
tion, a mass model has to be assumed [3] or the correlation 
with mass sensitive observables, most notably the muon 
number in the deep detector, has to be exploited for an es- 
sentially two-dimensional unfolding [4] . 

The energy resolution improves with energy and ap- 
proaches 0.05 in logio£^, or 12% in E, at about lOPeV 
for zenith angles less than 30°. The angular resolution is 
better than 1°, almost independent of energy and zenith an- 
gle. The detector coverage is ^ w 3 km^ sr for IceTop 
alone and Afl w 0.3 km^ sr for coincidences with the in- 
ice detector. 



4 First results 

Energy spectrum: The shower reconstruction from Ice- 
Top signals has been developed mainly using data taken 
with 26 stations (nearly 1/3 of the complete detector) in 
2007. The spectra of the shower size parameter 5i25 for 
three different zenith angular ranges are shown in Fig. 3. 
Except for the threshold region, this parameter is a close 
proxy for the primary energy for a given zenith angle range. 
The relation between 5i25 and the true energy is mass de- 
pendent. Under the assumption of an isotropic cosmic ray 
flux, the 5i25 spectra for different zenith angles should 
yield the same energy spectrum. It has been shown that 
this can only be achieved under the assumption of a mixed 
composition [3]. 

A first evaluation of IceTop data from the 2010 season with 
79 IceCube strings and 73 IceTop stations is reported else- 
where in these proceedings [5]. 

Mass composition using IceTop and deep-ice coinci- 
dences: As emphasized in the introduction, a strength of 
IceCube is the possibility to measure high energy muons in 
the deep ice in coincidence with the shower reconstructed 
in IceTop. The first analysis of such coincident data is pre- 
sented at this conference [4]. The data set is constrained to 
a small fraction of the detector and a relatively short time 
period (about 1 month). However, the results are currently 
more affected by systematic uncertainties than by statistics. 

Figure 4 shows a simulation of the correlation between the 
parameter Kro, which measures the muon energy in the 
deep ice, and the shower size parameter 5i25, which is a 
measure of the electromagnetic component of the shower. 
Different primary masses populate different bands in this 
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Figure 3: Measured spectrum of the shower size parameter 
5i25 for three different zenith angular ranges (data from 6 
months with 26 IceTop stations). 
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Figure 4: Mass composition measurement: plotted is the 
muon energy parameter K^Q measured in the deep ice ver- 
sus the shower size 5*125. The band is obtained by sim- 
ulating proton and iron primaries. The shading indicates 
the proton content, with 100% protons at the bottom und 
100% iron at the top. Indicated are also curves of constant 
primary energy as determined by a neural net (the labels are 
in units logio E/GeV). The points and the vertical bars are 
averages and dispersions, respectively, of the measured ex- 
perimental distributions of a logio ^1^5 ^™ (indicated by 
the horizontal bars) projected unto the logj^g Xyo axis. 




Figure 5: IceCube event display showing light signals in 
DOMs: Candidate muon bundle with a high-px muon (on 
the right). 
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Figure 6: Muon counting in IceTop: Distribution of tank 
signals for various cuts on the signal expectation Sexp in 
the energy range between 1 and 30 PeV. 

plot. Indicated are also curves of constant primary energy. 
The data points in log £'125 bins are plotted with the disper- 
sion of their distribution along the K^q axis. 

With a neural network, the primary energy and a parame- 
ter related to In A (A is the atomic mass number) has been 
determined yielding an energy and average In A spectrum 
between 1 and 50 PeV. The method allows for the extrac- 
tion of multiple mass contributions from fits to the neural 
network output, see details in ref [4]. 

High-px muons: Prompt decays of heavy flavour 
hadrons occurring in the first interactions are expected to 
produce muons with large transverse momentum. The pre- 
dictions are still very model dependent. In these proceed- 
ings, an analysis is presented [6] where high-px muons 
have been found as single tracks separated from a muon 
bundle by more than 200 m (Fig. 5). Current work concen- 
trates on understanding the systematic uncertainties in the 
resulting px distribution. 

Searching for PeV gamma rays: IceCube can effi- 
ciently distinguish PeV gamma rays from the background 
of cosmic rays by exploiting coincident in-ice signals as 
veto. Gamma-ray air showers have a much lower muon 
content than cosmic ray air showers of the same energy. 
Candidate events are selected from those showers that lack 
a signal from a muon bundle in the deep ice. Results of 
one year of data, taken in the 2008/2009 season when the 
detector consisted of 40 strings and 40 surface stations, are 
presented at this conference [7]. The projected gamma-ray 
sensitivity of the final detector is also given. 

Muon counting in IceTop: The muon content of a 
shower is a mass sensitive observable since heavier pri- 
maries tend to have a higher muon abundance. Although 
the number of high energy muons in the muon bundle near 
the shower core is most sensitive, muons counted at the 
surface (at typically much lower energies than in the deep 
ice) provide additional information on the mass. The com- 
parison of both methods allows one to test hadronisation 
models. 
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Figure 7: Average scaler rates of several tanks in IceTop 
during a period in February 201 1 when a Forbush decrease 
occured. 

Muonic and electromagnetic signals in IceTop can in gen- 
eral not be distinguished. However, muons show up as a 
relatively constant signal of about 1 VEM at larger dis- 
tances from the shower axis where the expectation value 
of a tank signal, Sexp, becomes small compared to a muon 
signal (Sexp is obtained from the fit to the lateral shower 
distribution). As shown in Fig. 6, the muon signal becomes 
more and more prominent when requiring smaller Sexp- 
For Sexp < 0.125 VEM the figure illustrates that the num- 
ber of muons can be well fitted. Comparing these muon 
numbers as a function of energy to simulations of different 
primary masses an independent information on the mass 
composition is obtained. 

Heliospheric physics: The IceTop tanks detect sec- 
ondary particles produced by cosmic rays in the multi-GeV 
energy regime interacting in the atmosphere with a count- 
ing rate exceeding 1 kHz per detector With IceTop, he- 
liosperic disturbances of this rate can be studied with very 
good time resolution. Since each detector has a different 
threshold setting, it is also possible to estimate the en- 
ergy spectrum of the cosmic rays related to such events. 
In an IceCube contribution to this conference [8] the per- 
formance during a Forbush decrease observed in February 
2011 is demonstrated (Fig. 7). 

Other cosmic ray results: At this conference the Ice- 
cube Collaboration also reports results on the observation 
of cosmic ray anisotropics in the southern sky [9]. These 
results and the cosmic ray studies in ref. [10] use all muons 
detected in the deep detector which gives a much larger an- 
gular coverage. 
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Abstract: The IceCube Observatory at the South Pole is composed of a deep detector and a surface detector, IceTop, 
both of which use Cherenkov Hght to detect charged particles. Cosmic ray air showers contain muhiple particle 
components: in particular, electrons and muons detectable at the surface by IceTop, and high-energy muons detectable 
by the deep IceCube detector, in relative amounts that depend on the primary cosmic ray mass. Thus, coincident events 
can be used to measure both the energy and the mass composition. Here, a neural network is trained with simulations to 
map observables from the two detectors (input) into energy and mass estimators (output). Experimental data is then run 
through the same network, to measure the energy spectrum and average logarithmic mass of cosmic rays in the energy 
range of about 1-30 PeV. 
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1 Introduction 



2 Reconstruction 



Measuring more than one particle component of a cosmic 
ray air shower is a powerful tool for separating light and 
heavy nuclei in Extensive Air Shower (EAS) data, at ener- 
gies at the knee and above where direct measurements of 
cosmic rays are not possible. The IceCube Observatory at 
the South Pole is being used in such a manner, with an ar- 
ray of Ught sensors (Digital Optical Modules, or DOMs) 
buried on strings between 1450 and 2450 m (herein "Ice- 
Cube") together in coincidence with a corresponding array 
of DOMs in frozen water tanks on the surface ("IceTop"). 
An IceTop "station" is two tanks separated by 10 m; Ice- 
Cube strings and IceTop stations are separated by 125 m 
[1]. 

Both instruments measure the Cherenkov light emitted by 
charged particles through the surrounding medium. In Ice- 
Top, the medium is water frozen in tanks at the site. The 
DOMs in the tanks measure light from the electromag- 
netic and GeV muonic components of the EAS. In IceCube, 
the medium is the Antarctic icecap, and the deeply-buried 
DOMs measure light from high-energy (TeV) muons bun- 
dled near the central axis of the shower. 



Events in IceTop are reconstructed by a likelihood method 

[3], comparing the detected signal locations, charges, and 
times from hit stations (as well as the locations of not-hit 
stations) to what is expected from a cosmic ray EAS. Sig- 
nal times are compared to an expected timing profile, and 
signal charges are compared to an expected lateral distri- 
bution function (LDF). The LDF used, a function of the 
perpendicular distance from the shower axis, r, is known 
as the "Double Logarithmic Parabola,": 



S{r) = 5ref 



( — 



-/3-'-logio(:B^) 



(1) 



The logarithm of the signals S are assumed to have nor- 
mal distributions, and are expressed in units of "vertical 
equivalent muons" (VEM). Here, i?ref represents a "refer- 
ence distance" from the shower axis, and S^d is the signal 
strength at that reference distance. The reference distance 
found to have the most robust measurement at these ener- 
gies is 125 m [3], so 5125 is the observable representing the 
"shower size". 

Events in IceCube are reconstructed with a similar philos- 
ophy: the signals NpE, measured in photoelectrons, are 
compared to an LDF which is a function of perpendicular 
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distance to the track d (described in detail in [4, 5, 6]): 

^J{cic^{z) ■ Xo)d 
(2) 

This function of d is dominated by the decaying exponen- 
tial; the "slope" of this exponential is the attenuation length 
of light in the ice. Since the clarity of the ice changes due 
to well-measured horizontal dust layers |7, 8], this slope 
is treated as a bulk attenuation length Aq multiplied by a 
depth-dependent correction factor Cice based on scattering 
length data measured at different depths z. The first term, 
containing the slant depth from the snow surface X, cor- 
rects for the ranging-out of muons as they penetrate deeper 
into the ice; the parameters a, b, and -y^ in the first term are 
constants. 

The overall normalization of this LDF scales with the 
energy deposited by muons in the detector, and this is 
parametrized with K^q, which is the expected signal (given 
by Equation (2)) evaluated using reference values of X 
= 1950 m, Cice = 1> and at a perpendicular distance d = 
70 m. Since the track direction affects the signal expec- 
tation NpE at all the DOMs (by changing d and X), the 
same UkeUhood function can also be used to find the track 
direction. These two reconstruction techniques (using sur- 
face signals to find the core position and S'125 in IceTop, 
and using signals in IceCube to find the track direction and 
K^o) are used iteratively to find a reUable best-fit track. 

3 Data, Simulation, and Event Selection 

This analysis used data from August 2008, when the de- 
tector was in its 40-string/40-station configuration, for 
an overall detector livetime of 29.78 days. Because of 
imsimulated effects near the detector threshold due to 
snow buildup over tanks deployed before 2007, data from 
the subarray of IceTop deployed after 2007 were used. 
Monte Carlo simulated events were produced using the 
CORSIKA air shower generator [9] with the SIBYLL- 
2.1/FLUKA-2008 hadronic interaction models [10, 11], 
and an atmospheric model representing austral winter at 
the South Pole. Five particle species (proton, heUum, oxy- 
gen, silicon, and iron) were generated according to an 
spectrum from 1 TeV to 46.4 PeV. The showers (3000 of 
them per species per third of a decade in energy) were gen- 
erated uniformly over all azimuths and to a zenith angle 
of 65°, oversampled 100 times, and thrown over a circle 
of radius 1200 m centered on IceTop. The response of the 
IceTop tanks, the propagation of the high-energy muons 
through the ice to the depth of IceCube, the Cherenkov 
photons propagating to the DOMs, and the response of the 
DOMs themselves, are included in a detector simulation. 

QuaUty cuts were applied to select those events which were 
well-reconstructed and contained within both the IceTop 
area and IceCube detector volume. Parameters such as 
track length and reconstructed effective propagation length 
(Ao) were used to quantify reconstruction quality (see [6] 



for details). These cuts allow for a resolution (containing 
68% of events) of better than 9 m in core position and less 
than 0.5° in track direction. The final event sample contains 
239797 events from the August 2008 experimental data and 
20289 total simulated events of five primary species. 

4 Neural Network Mapping Technique 

Once S'125 and K-^o have been reconstructed and quaUty 
cuts applied, simulations show that light and heavy nu- 
clei are separated in this two-dimensional parameter space, 
as shown in Figure 1. The relation between the K^q- 
5*125 space and the mass-energy space is non-linear, there- 
fore a mapping technique is required to correlate one space 
to the other. 




Figure 1: Fraction of each bin populated by protons 
(when only protons and iron are included in the sample). 
Dark grey indicates 100% iron, Ught grey indicates 100% 
protons; intermediate greys indicate overlapping popula- 
tions. The dotted black hues approximating energy con- 
tours guide the eye for both nuclei. 

A neural network was chosen for this work, which consists 
of a set of input parameters (in this case, K-^o and ^125) 
which are connected to a set of output parameters (in this 
case, logio(£') and ln{A)) through a series of nodes which 
are arranged in layers. Each node is connected to other 
nodes in the previous and subsequent layer via a series of 
weights. At each node, an activation function acts on its 
input parameters as modified by the weights. Both the two 
inputs and the two outputs are renormalized so that they are 
numbers between zero and one. 

The weights relating the inputs to the outputs are deter- 
mined by "training" the network on a subsample of the 
Monte Carlo simulations (1/4 of the events) for which the 
true energy and true mass are known. Through a series of 
learning cycles, the network adjusts the weights to improve 
the accuracy of reconstructed outputs. A "testing" sample 
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(an independent 1/4 of the events) are put through the net- 
work after each learning cycle and the errors monitored to 
ensure that the network is not becoming too specific to the 
training events, or being "overtrained". The remaining half 
of the simulated data is known as the "analysis" sample and 
was used for the final steps described in section 6 below. 

5 Energy Spectrum 

For each event of experimental data, the neural network as- 
signs a reconstructed energy and a reconstructed mass pa- 
rameter Reconstructed energies agree well with the true 
energies of simulated events; Figure 2 shows the energy 
resolution, and the reconstruction bias, as a function of en- 
ergy. 
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where 77 is the efficiency (the ratio of simulated events left 
after all cuts to the number generated, which is a function 
of £'0), A is the area over which the CORSIKA showers 
were thrown, r is the livetime of the detector (in this case 
29.78 days), and fl is the solid angle over which the events 
were generated. The resulting energy spectrum is shown in 
Figure 3. 
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Figure 3: Measured all-particle flux as a function of particle 
energy, from August 2008 experimental data, with statisti- 
cal error bars. Triangles indicate alternate spectra derived 
from data with K70's shifted up or down by 22.2% (see 
section 7). 
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Figure 2: Upper: The energy resolution of the neural 
network output (sigma of distribution of log{Ereco) — 
log{Etrue))- Lower: The bias or misreconstraction of en- 
ergy (mean of that distribution). 

From these reconstructed energies, one can create an all- 
particle energy spectrum from experimental data. For a 
given flux ^1/ as a function of energy Eq, 



1 dN 



1 0.4343 dN 



riAnr dE riAnr Eq cflogio(-Eo)' 



(3) 



6 Mass Composition 

Within each slice in energy and for each simulated species, 
the neural network produces a distribution of mass outputs 
which is called a "template histogram". Examples of tem- 
plate histograms for three kinds of primaries are shown in 
Figure 4, for one slice in energy as an example. Data, when 
put through the same neural network, also has a histogram 
of outputs which can be decomposed into a linear combi- 
nation of the template histograms of the individual species 
(proton, oxygen, etc.). A minimizer finds the optimal mix- 
ture of simulated species to match the data. 

Proton and iron template histograms alone are not suffi- 
cient to reproduce the data, and so intermediate nuclei are 
necessary. Because there is a great deal of overlap between 
the template histograms of all five nuclei, we matched the 
data to a combination of three template histograms: pro- 
tons, iron, and "intermediate nuclei" (which is a 50-50 mix- 
ture of helium and oxygen). When this procedure is applied 
to a "hand-mixed" sample of Monte Carlo events (treating 
the sample like data, with different \n{A) at different ener- 
gies), the fit found mixtures which correctly reproduced the 
In (A) at all energies. When the procedure was applied to 
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Helium+Oxygen 

Iron 

Data 




NN Output: Type 

Figure 4: Three "template histograms": protons, iron, and 
a 50-50 mixture of helium and oxygen. The energy bin 
shown here is logiQ(£'/GeV) from 6.9 to 7.1. 

experimental data from August 2008, the In(^) as a func- 
tion of energy was computed and is shown in Figure 5. 



Preliminary 



• Ice Model: Since K-;q depends on the propagation 
of photons through the ice, the effect of two differ- 
ent models of the attenuation of light in dust layers 
(known as AHA [7] and SPICE [8]) were investi- 
gated. Percent error in Kiq: 9. 7% 

• DOM Efficiency: The efficiency of the DOMs in Ice- 
Cube depends upon a number of factors, which have 
been measured in a controlled setting [14]. Percent 
error in Kiq: 8% 

• Errors introduced from corrections: The Monte 
Carlo simulation used for this work did not include 
two relevant effects: a recent improvement of the pa- 
rameterization of the light yield from muon bundles 
in IceCube, and the accumulation of snow on the sur- 
face above IceTop. The effects were quantified and 
corrected for, but some error due to the application of 
this correction was estimated. Percent error in K-jq: 
14.8% and 4.7%, respectively 

The total systematic shift in Kjq (added in quadrature) is 
22.2%. The effect of these systematics on the final results 
(i.e. the spectrum in Figure 3 and the mass composition in 
Figure 5) were studied by applying these errors as shifts to 
the input parameters of the neural network (i.e. K-jq) before 
it goes through the neural network, and seeing how the out- 
put parameters, and any analysis downstream of the neural 
network, changes. Both of these figures show alternative 
results with such a shift (both up and down) applied. If 
improved simulations eliminate the errors introduced from 
corrections, the total systematic error could be substantially 
reduced. 
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Figure 5: Measured In(^) as a function of reconstructed 
energy, with statistical error bars. Triangles indicate alter- 
nate results if K70 is shifted up or down by 22.2%. 



7 Systematic Errors 

Both Figures 3 and 5 include estimates of the systematic 
error due to a variety of effects: 

• Hadronic Interaction Model: Samples of events sim- 
ulated with both EPOS-1.99 [12] and QGSJET-II-03 
[13] were generated for comparison with SIBYLL. 
Percent error in JCto-' 9.7% 
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Abstract: The high statistics of cosmic ray induced muon events detected by the IceCube Observatory makes it possible 
to study the correlation of muon intensity with the stratospheric temperature over Antarctica with high precision. Using 
150 billion events collected by IceCube experiment over 4 years, the muon rate was found to be highly correlated with 
daily variations of the stratospheric temperature and exhibits a ± 8% aimual modulation. The correlation between the 
muon rate and the upper atmospheric temperature is related to the relative contribution of tt and K to the production of 
muons. Therefore it is possible to estimate the K/tt ratio from the seasonal variation of the muon rate, which was found 
to be 0.09 ± 0.04 at cosmic ray median energy of about 20 TeV. 
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1 Introduction 

When cosmic ray particles enter the Earth's atmosphere, 
they generate a hadronic cascade in which mesons are pro- 
duced, primarily pions and kaons. These mesons can either 
interact again or decay into muons. The relative probability 
of decay or interaction depends on the local density of the 
atmosphere, which in turn depends on the temperature [1]. 
The differential flux of muons with energies larger than 100 
GeV can be described with good approximation as [2] 



where (/)jv(-E^) is the primary spectrum of nucleons (N) 
evaluated at the energy of the muon. The first term in 
Eq. 1 corresponds to muon production from leptonic and 
semileptonic decays of pions, while the second term is re- 
lated to kaons. The constants and Axfi include the 
branching ratio for meson decay into muons, the spectrum 
weighted moments of the cross section for a nucleon to pro- 
duce secondary mesons, and those of the meson decay dis- 
tribution. The denominators in Eq. 1 reflects the competi- 
tion between decay and interaction of secondary mesons in 
the atmosphere. When Et^^k < e^r./f/cos^*, the meson de- 
cay is the dominant process, and muons are produced with 
the same spectral index as the parent cosmic rays. At high 



energy meson interaction dominates and the corresponding 
muon spectrum becomes one power steeper than the pri- 
mary spectrum. The characteristic critical energies x 
at a given atmospheric depth are inversely proportional to 
the atmosphere's density at that point, and therefore are af- 
fected by temperature variations. In an isothermal approxi- 
mation of the atmosphere, the density is described by an ex- 
ponential with a scale height of ho w 6. 19 km (over Antarc- 
tica). The numerical value applies to the lower strato- 
sphere, where most of the muons are generated. In this 
approximation x are proportional to the atmosphere's 
temperature in the perfect gas state limit. At a mean at- 
mospheric temperature of To = 21 1°K the critical energies 
are 6.^ — 111 GeV and €k = 823 GeV. The dependence of 
the critical energies on temperature is the main source of 
the seasonal variation in muon rate. This modulation was 
studied by underground experiments such as MACRO [3], 
LVD [4] and MINOS [5], and by AMANDA [6] and Ice- 
Cube [7]. Here we update the analysis with four years of 
IceCube data with an emphasis on the systematic effects 
that can be studied with a very large amount of data. 

2 Temperature Correlation 

The relation between the variation of temperature and the 
variation of muon intensity at a particular energy and zenith 
angle can be expressed in terms of a theoretical correlation 9 
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coefficient calculated from Eq. 1 as 



dT 



(2) 



which depends explicitly on the characteristic critical ener- 
gies e^^^K- 

Measured rates depend on the convolution of the muon dif- 
ferential spectrum with the detector response, which de- 
pends on muon energy and zenith angle. To compare 
with measurements, it is therefore necessary to calculate 
a weighted correlation coefficient as 



JdE,^^iE^,9)A,ffiE^,0) 



(3) 



where Aeff{E^, 0) is the effective detector area obtained 
from simulation. Eq. 3 defines the correlation coefficient 
for a particular zenith angle 9. The total correlation co- 
efficient is then obtained by averaging a*^{0) over with 
a weight given by the observed event angular distribution. 
With this definition the variation in muon intensity is 
given by 



th 



AT, 



(4) 



where Tg/ / is the effective atmospheric temperature as de- 
fined below. Since the rate i?^ of observed muons is pro- 
portional to the incident muon intensity 7^, it is correlated 
with the effective temperature as well 



AT, 



{Teff) 



(5) 



where is the experimentally determined correlation 
coefficient. 

Since muon production occurs over an extended portion of 
the upper atmosphere and the temperature depends on al- 
titude, it is necessary to define a parameter referred to as 
effective temperature, in order to quantify the relationship 
between variations in temperature and those in measured 
muon rate. 

The pion and kaon terms in Eq. 1 are derived from the in- 
tegral over the atmospheric slant depth X (in glcw?) of the 
muon production spectrum P^{E^^ 9, X), which in turn is 
given by the probability distribution for meson decay to 
muons integrated over the parent meson spectrum [2]. 

The effective temperature as a function of muon energy 
and zenith angle is defined as the actual temperature profile 
weighted by the muon production spectrum 



J dXP^iE^,9,X)T{X) 
JdXP^{E^,9,X) 



where P^{E^, 9, X) is the sum of muon production spec- 
trum from pion and kaon contributions. The rationale for 
this definition is that the depth dependence of the muon 
production spectrum weights the temperature with the re- 
gions of the atmosphere where the meson decay to muons 
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Figure 1 : The differential weighting function used in the 
calculation of the Teff, as a function of atmospheric depth 
X (continuous line), superimposed on the average atmo- 
spheric temperature profile (dashed line). 



occurs. The critical energies that appear in the production 
spectrum are evaluated at the mean temperature. To, but 
the production profile does not depend strongly on temper- 
ature. 

This formulation of the effective temperature differs from 
that in Grashorn et al. [8], where the low energy limit of 
the temperature derivative of muon production spectrum 
produces an unphysical discontinuity. The IceCube obser- 
vatory is located at a depth of >1.3 km.w.e./cos6' so only 
muons with energies above about 400 GeV/cosfl can reach 
and trigger the detector. At these energies muon decay and 
energy loss in the atmosphere are negligible. This enables 
us to use analytic forms for muon production analogous 
to Eq. 1 without accounting for decay and energy loss of 
muons in the atmosphere and hence to obtain a physically 
correct result without any discontinuity. 

To compare predictions with measurements it is necessary 
to determine the convolution with the detector response 
function 



Teffie) = 



JE, ^dXP^{E,,9,X)A,ff{E^,9)T{X) 
fE^fdXP^{E^,9,X)A,jf{E^,9) ' 

(7) 

where the denominator, as in Eq. 3, is the total measured 
muon intensity. As for the correlation coefficient, the total 
effective temperature Te/ / is the weighted average of Eq. 7 
over the event zenith distribution. Fig. 1 shows the dif- 
ferential weighting function used in the calculation of the 
effective temperature along with the average atmospheric 
temperature profile over the seasons. It peaks at about 100 
g/crn^, which is where most of the muons are produced. 



(6) 3 Muon and Temperature Data 



In this analysis we used 150 billion muon events col- 
lected by the partially completed IceCube Observatory 
from March 2007 to April 2011 (see [9] for an overview 
on IceCube). These events are generated by cosmic rays 
with median energy of about 20 TeV. 
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Figure 2: The daily atmospheric temperature profiles over Antarctica produced by NASA AIRS instrument on board the 
Aqua satellite [10] are shown from 2007 to 201 1 (in color code), along with the relative modulation in the measured muon 
event rate (black line) and the relative variation of the effective temperature (grey line). Note that the statistical uncertain- 
ties in the data are between about 0.1 Hz (in 2007) and about 0.16 Hz (in 2011) on the daily rates. The corresponding 
uncertainties on the relative rate variations are 0.02 % and 0.008 %, respectively. 



The muon rate increased substantially over these four years 
as new detectors were added during each construction sea- 
son. As the instrumented volume increased, the probability 
that one data record included two or more separate cosmic- 
ray events increased from about 1% to about 4%. A cor- 
rection to the daily recorded rate was therefore applied to 
obtain a corrected rate of muon events, i?^. 

The atmospheric temperature profile data used in this anal- 
ysis were collected by the NASA Atmospheric Infrared 
Sounder (AIRS) on board the Aqua satellite. Daily atmo- 
spheric temperatures at 20 different pressure levels from 1 
to 600hPa above the South Pole were obtained from the 
AIRS Level 3 Daily Gridded Product available on NASA 
Goddard Earth Sciences, Data and Information Services 
Center (GES DISC) [10]. Using these data the daily effec- 
tive temperature T^ff was calculated based on the zenith- 
weighted average of Eq. 7. 

4 Results and Determination of K/tt Ratio 

Fig. 2 shows the measured as a black continuous line 

AT f f 

along with ^j, as a black dashed line. The figure also 
shows the actual atmospheric temperature profile as a func- 
tion of pressure level (equivalent to atmospheric depth X). 
The statistical uncertainties in the measured muon rate are 
too small to show in the figure. Note that besides the large 
seasonal modulation, the daily muon rate is strongly cor- 
related with short time temperature variations in the upper 
atmosphere. 

Based on Eq. 5, the experimental temperature coefficient 
was determined from regression analysis and found to be 
a^^P = 0.860 ± 0.002 (stat.) ± 0.010 (syst.). The experi- 
mental systematic uncertainty on a'^^ is dominated by the 
effective area A^f / (S^ , 0), which is used in the calculation 
of the effective temperature and of the theoretical correla- 
tion coefficient in Eq. 3. Most of the detected muons range 
out within the large instrumented volume of IceCube, and 



the energy profile of the effective area depends on the dis- 
tribution of depths the muons reach within the array. The 
spread on this distribution, translates into an estimated un- 
certainty in the experimental correlation coefficient of 0.01 . 

Since the temperature correlation coefficient depends on 
the relative contribution of pions and kaons, it is possible 
to use the seasonal variations of the muon rate to determine 
the K/tt ratio. 

The effective temperature Tg/ / is relatively insensitive to 
variations in the cosmic ray spectral index, the proton at- 
tenuation length, the critical energies and K/tt ratio because 
the dependence cancels to a large extent due to the normal- 
ization in Eq. 7. 

The theoretical correlation coefficient a^, on the other 
hand, depends primarily on the critical energies and on the 
K/tt ratio. Changing the cosmic ray spectral index and pro- 
ton attenuation length within a wide range, has an effect 
smaller than 1%. Therefore it is possible to use the param- 
eters for attenuation lengths and spectrum weighted mo- 
ments from Ref. [2], assuming a cosmic ray spectral index 
of -2.7. The critical energies evaluated at the average effec- 
tive temperature of To = 21 1°K are used. 

In particular, the kaon to pion ratio Rk/tt — de- 
pends on the spectrum weighted moments Z^k and Z^r^ 
of the cross section for a nucleon N to produce secondary 
kaons and pions, respectively, from a target nucleus in the 
atmosphere. The dependence on the spectrum weighted 
moments Znk and Z^r^ is implicit in the parameters 
and AKf_i in Eqs. 1, 3, 7. 

The nominal value of K/tt ratio is taken to be [2] 

= =0.149 ±0.060, (8) 

which is based on laboratory measurements below 100 
GeV center of mass energy. The 40% uncertainty corre- 
sponds to that in the current cosmic ray interaction mod- 
els [11]. 
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Figure 3: The preliminary determined by IceCube as 
a function of the K/tt ratio (continuous Une), with its error 
given by the band, and the predicted (dashed Hne) with 
the corresponding error given by the band. The intersection 
region between the two bands is at = 0.09 ± 0.04. 

By calculating the theoretical correlation coefficient as 
a function of Rk/-k ^nd comparing it with the experimental 
value, it is possible to measure the kaon to pion ratio for 
proton interaction with atmospheric nuclei (mainly nitro- 
gen) at cosmic ray particle energy of 20 TeV. 

Fig. 3 shows the experimental and theoretical values of ax 
as a function of K/tt ratio for the IceCube data. While, as 
mentioned above, a!^^ is almost insensitive to RK/-n, the 
value of a^rp has a strong dependence. Since the statistical 
uncertainties are very small, the band on is dominated 
by systematic uncertainties. The band on atj^ reflects the 
theoretical uncertainties in the calculation of the correlation 
coefficient mentioned above, and conservatively estimated 
to be 1%. The crossover in Fig. 3 is at Rk/^, = 0.09 ± 0.04. 
Fig. 4 shows a comparison of the IceCube measurement 
with other observations. 
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Figure 4: A compilation of selected measurements of K/tt 
for various center of mass energies. Data points are from 
NA49 [12, 13],E735 [14], STAR [15] and MINOS [5]. The 
horizontal line and gray band represents the reference value 
K/tt ratio 0.149 ± 0.060 [2, 1 1]. 

the atmospheric muon charge ratio from MINOS [16] and 
OPERA [17]. 

6 Acknowledgements 

The temperature data used in this study were acquired 
as part of the NASA's Earth-Sun System Division and 
archived and distributed by the Goddard Earth Sciences 
(GES) Data and Information Services Center (DISC) Dis- 
tributed Active Archive Center (DAAC). 

References 



5 Conclusions 

Using 150 billion cosmic ray induced muon events col- 
lected in four years by IceCube, a strong correlation of 
the daily observed muon rate with the stratospheric tem- 
perature was observed, along with a ±8% annual modu- 
lation. The K/tt ratio at 20 TeV cosmic ray energies was 
determined by comparing the observed temperature corre- 
lation coefficient with the theoretical one, and found to be 
Rk/^ = 0.09 ± 0.04. 

The value obtained with IceCube implies that Z^k ^ 
0.0071, which is about 40% lower than its nominal value 
from Ref. [2]. In calculating the theoretical correlation co- 
efficient the sum Zjssk + Zntt was kept constant to its nom- 
inal value 0.0908. One way to reconcile the measurement 
of Rk/h with other results is by reducing the amount of 
associated production pA^ nK^ . Keeping the nominal 
value of Z^K- — 0.0028, while reducing Zpfc+ from its 
nominal value of 0.0090 to 0.0043 is also likely to give a 
better agreement with recent measurements of the ratio of 



[I] RH. Barret, et ai. Rev. Mod. Phys., 1952, 24(3), 133. 
[2] T.K. Gaisser, 1990, Cosmic Rays and Particle Physics, 

ed. CUP 

[3] M. Ambrosio, et ai, Astrop. Phys., 1997, 7, 109. 
[4] M. Selvi, Proc. of 3P* ICRC, 2009, Lodz, Poland. 
[5] P Adamson, et ai, Phys.Rev., 2010, D81,01200. 
[6] A. Bouchta, Proc. of 26*'* ICRC, 1999, Salt Lake City, 
UT, U.S.A. 

[7] S. Tilav, etai, Proc. of 3P* ICRC, 2009, Lodz, Poland, 
arXiv: 1001. 0776. 

[8] E.W. Grashom, et ai, Astrop. Phys., 2010, 33, 140. 

[9] H. Kolanoski, IceCube summary talk, these proceed- 
ings. 

[10] http://disc.sci.gsfc.nasa.gov/AIRS/data-holdings. 

[II] G.D. Barr, et ai, Phys. Rev., 2006, D74, 094009. 
[12] S.V. Afanasiev, et ai, Phys. Rev., 2004, C66, 054902. 
[13] C. Alt, et ai, Em. Phys. J. C, 2006, 45, 343. 

[14] T. Alexopoulos, et ai, Phys. Rev., 1993, D48, 984. 
[15] C. Adler, et al, Phys .Lett., 2004, B595, 143. 
[16] P Adamson, et ai, Phys. Rev., 2007, D76, 052003. 
[17] N.Agafonova,efflZ., 2010, ai-Xiv:1003. 1907. 12 



32nd International Cosmic Ray Conference, Beijing 201 1 




Measurements of the Air Shower Parameters with IceTop 

The IceCube Collaboration^ 

^See special sections in these proceedings 



Abstract: We study the lateral distribution function (LDF) of signals in the IceTop air shower detector as a function 
of distance from the air shower core. The completed IceTop detector consists of 81 stations with two tanks each. It can 
now study the signals at distances approaching 1 km from the core position. We discuss the general shape of lateral 
distributions of the signal and its dependence on the shower zenith angle and primary mass. We also show the simulated 
individual tank signal lateral distribution for a large number of simulated proton and iron showers. We find that the form 
of the lateral distribution function used for more widely spaced arrays of water Cherenkov detectors, Haverah Park and 
Auger in the EeV range, can also be used with appropriate parameters to describe IceTop data in the 10-100 PeV range. 

Corresponding authors: Shahid Hussain^ (shahid@bartol.udel.edu), Todor Stanev^ {stanev@bartol.udel.edu). Scrap 
Tilav^ (tilav@udel.edu) 

'^University of Delaware, Newark, DE 19716, U.S.A. 
Keywords: IceCube, IceTop, lateral distributions 



1 Introduction 

IceTop, the surface air shower array above the IceCube 
neutrino detector, was completed in the 2010/2011 deploy- 
ment season. IceTop consists of 81 stations, each of which 
has 2 tanks of area 2.7 m^ containing 90 cm of clear ice. 
The average distance between the tanks in a station is 10 m 
and the average distance between stations is 125 m. The 
Cherenkov light generated by the charged particles that 
hit the tanks is collected by two digital optical modules 
(DOMs) that run at different gains to increase the dynamic 
range. The signal strength is measured in vertical equiva- 
lent muon units (VEM), i.e. the signal that a 1 GeV vertical 
muon produces in the tank. A station triggers when there 
are signals above the threshold (0.16 VEM) in both tanks 
within one microsecond. In this way we avoid triggering 
on coincidental muons that belong to different atmospheric 
cascades. 

The current shower reconstruction of IceTop events is 
based on a procedure which was designed when the array 
contained only 26 stations in a much smaller area [1]; the 
reconstruction procedure is applied to showers that trigger 
at least five stations. The signal lateral distribution func- 
tion in this procedure does not include tanks with zero sig- 
nal (i.e. tanks that do not trigger) and the fitting routine 
accounts for these tanks with a separate, no-hit probability, 
term in the likelihood function. 



In this paper we instead use a lateral distribution function 
of the form [2] 

Sir) = Ax r-^n+r/ro)^ (1) 

where S{r) is the signal at a perpendicular distance r 
from the shower core in shower coordinates; r/, rp, and A 
are fit parameters. This form has been used for Haverah 
Park [3] and Auger [4J to fit showers observed with wa- 
ter Cherenkov tanks. IceTop is at a much higher altitude 
(2835 m) than HP (sea level) and significantly higher than 
Auger (1400 m) and collects data in a different primary en- 
ergy range with detectors spaced by 125 m as compared to 
several hundred meters for HP and 1500 m for Auger. We 
investigate here the extent to which the lateral distribution 
form used for the other Cherenkov shower detectors scale 
to the location and energy range of IceTop. Everywhere in 
this paper, S(r) gives the signal strength at a perpendicu- 
lar distance r from the shower core in shower coordinates. 
We include tanks with and without a signal directly in the 
lateral fit, both for simulated and observed air showers. 

2 Monte Carlo Calculation 

To study the lateral distribution we have simulated air 
showers initiated by protons and iron nuclei with fixed pri- 
mary energies of 10 and 100 PeV and fixed zenith angles 
of 0, 25, and 45 degrees. Air showers are simulated with 
CORSIKA-SIBYLL [5] and the detector simulation uses 
Geant4 [6] for the tank response. 50 showers per primary 
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Figure 1: LDF of the signal strengths (weighted by square 
of the distance from core) for 100 PeV simulated proton 
and iron showers at three zenith angles. 
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Figure 2: Ratio of the signal strength in simulated iron to 
proton showers for three zenith angles. 

type, energy, and zenith angle were dropped on the IceTop 
array 100 times each within a 600 meter radius from the ar- 
ray center. As a result, there are 5,000 showers in each set. 
The simulated showers were reconstructed with the current 
standard procedure and all results presented below use the 
reconstructed shower core position, direction, and shower 
energy for event selection. Fewer than 200 showers in each 
set were not reconstructed well and are not analyzed here. 

Figure 1 shows the average LDF of the IceTop signals for 
simulated proton and iron initiated showers at the three 
zenith angles. As expected, the proton showers have higher 
signals close to the shower core and iron showers have 
higher signal density at large distances from the core^ 
The intersection point of the proton and iron signals LDF 
changes with zenith angle and increases significantly even 
at the modest zenith angle of 45 degrees. The ratio of iron 
to proton signal LDF is presented in Figure 2. While for 
strictly vertical showers the intersection point is at about 
50 m from the shower core; it is between 120 and 150 m 



Figure 3: The ratio of RMS of the signal to its strength in 
simulated 100 PeV proton showers at three different zenith 
angles. 

away from the core for ^=25 degrees and above 300 m at 
45 degrees. Furthermore, although not shown here, the in- 
tersection point for 10 PeV showers is relatively larger for 
each angle; it is around 130, 180, and 400 m for 6»=0, 25, 
and 45 degrees, respectively. The IceTop reconstruction 
procedure currently uses the signal at 125 m, 5*125, from the 
shower core as the energy-related parameter, which seems 
to be mass independent for showers near the peak of the 
angular response for IceTop (25 degrees) as the classical 
papers of A.M. Hillas [7] recommend. However, as we see 
in Figure 2 for IceTop, there is no single distance from the 
shower core that is independent of mass for all energies 
and angles. Minimizing the fluctuations is also desirable 
in the choice of an optimum distance r used for energy as- 
signment. Figure 3 shows the ratio of RMSs/S'(r) for E = 
100 PeV proton showers and the three angles. The signal 
fluctuations are almost constant between 50 m and 150 m 
from the shower core for all three angles. The fluctuations 
are higher for the most inclined showers and still constant 
below 150 m distance from the core, the goal is to apply 
the lateral distribution function to experimental data which 
may include a mixture of proton, iron, and several other 
primaries. Therefore, it is instructive to look at the VEM 
range of signals for a given primary energy for both proton 
and iron together As an example, we show in Fig. 4 the 
signals measured by IceTop in 100 proton and 100 simu- 
lated showers of primary iron nuclei. At distances around 
100 m the signal variation is sHghtly larger than a factor of 
two. At much larger distances, around 500 m, the signal 
variation increases to more than one order of magnitude. 
A fraction of these fluctuations is due to the fact that we 
plot proton and iron induced showers together. There is a 
strong increase in fluctuations for each species at distances 
greater than 300 m, as shown in Fig. 3. Note also the rela- 
tively small statistics at distances below 100 m. 

1. Everywhere in this paper, distance from the core is the per- 
pendicular distance from the core in shower coordinates. 
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1280 and 1474 m. While rj values are similar for proton and 
iron showers, the ro values are always significantly smaller 
by about 200 m than those of iron showers. The respective 
normahzation parameter A is 8.58 x 10^ , 7.29 x 10^, and 
2.02 X lO*" VEM for proton showers. 

One has to note that the fits are not very good for core dis- 
tances less than about 50 m. The reason is that we are using 
reconstructed position of the shower core. The average er- 
ror in the core position is less than 20 m but this error still 
affects the LDF at small distances in a negative way. The 
statistics at small distances is also low. For these reasons 
the measured or interpolated signal strengths at small dis- 
tances are not reliable. 



Figure 4: The signals recorded by all non-zero tanks for 
100 proton and 100 iron induced vertical showers of energy 
100 PeV are plotted as a function of the distance to shower 
core. 
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Figure 5: The average LDF of the signal strength for 100 
PeV simulated iron showers at three zenith angles fitted 
with the HP-like function. 

3 Fitting the simulated signal LDF 

Both HP and Auger have fitted their data with a lateral dis- 
tribution function similar to Eq. 1 ro fixed. Here we fit the 
average of the simulated showers at each energy with the 
same form. As an illustration, the signal LDF of the simu- 
lated 100 PeV iron showers and the fits obtained are shown 
in Fig. 5. As expected, the lateral distribution becomes 
flatter with angle. The parameter rj is 2.22, 2.21, and 2.06 
for 0, 25, and 45 degrees, respectively. The respective ro 
values are 1410, 1507 and 1783 m. The respective normal- 
ization parameter A is 7.18x10^, 5.31 xIQS, and 1.25x10^ 
VEM. The values of all fits are smaller than 1 per degree 
of freedom. 

The fits of simulated proton showers have similar results. 

The T] values are 2.25, 2.22, and 2.10 for 0, 25, and 45 
degrees, respectively. The respective ro values are 1215, 



4 Fitting of individual showers 

The next step in the study of the signal lateral distribution 
function is the fitting of individual showers. We have at- 
tempted to fit individual simulated and experimentally de- 
tected showers with the function of Eq. 1 . The fitting pro- 
cedure is less stable when applied to individual showers 
especially because there are usually few points at distances 
smaller than 100 m from the shower core. This makes the 
ro parameter vary even more than in the case of average 
lateral distribution from a large number of showers. 

The well fitted showers, however, show a good agreement 
between the experimentally detected and the simulated 
showers. To compare these two sets of showers we chose 
experimental showers with standard IceTop reconstructed 
parameters very close to the simulated (fixed primary en- 
ergy Ep and zenith angle 6) ones. For vertical showers, for 
example, we chose showers with cos(^) > 0.95 and with 
1.97 < logio(£'p/PcV) < 2.03. Figure 6 shows the LDF 
fit of a simulated proton (top) and an experimentally de- 
tected shower (bottom). Both showers have a large number 
of triggering stations - 48 stations in the simulated shower 
and 46 in the experimental one. The experimental shower 
shown in Fig. 6 is from 2010 when IceTop consisted of 73 
stations (or 146 tanks). The tanks with '0' signal are not 
shown, but they are included in the average for each annu- 
lar bin of r. In the graph of the experimental shower we 
also show the average signal strength calculated in loga- 
rithmic radial bins around the reconstructed shower core. 

The fits have a of 0.6 p.d.f. However, the performance 
of the fit will depend on x^ distributions from fits on a rea- 
sonable sample of showers; a detailed study is needed in 
this regard. 

5 Summary 

We have studied the lateral distribution of air shower sig- 
nals in IceTop tanks using a function similar to the one used 
by Haverah Park and Auger Observatory. Tanks with sig- 
nals below threshold are included as zeros in averaging the 
signal in each radial bin. The simulated proton and iron 
showers at 10 and 100 PeV can be described with lateral 
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energy estimation reference distance, with the energy and 
zenith angle of the shower using an iterative procedure. In 
this regard, it is encouraging that, for showers in the peak 
of the angular distribution for IceTop (0 ^25 degrees) and 
energies 10-100 PeV, the variation of the iron to proton sig- 
nal ratio is around 1 and it has about 10% variation for the 
reference distance range of 100-150 m. It is also encour- 
aging that the signal fluctuations are almost constant in the 
distance range 50-150 m. 
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Figure 6: Example fits on two individual showers; both 
showers are almost vertical (cos(6')>0.95) with energies 
close to 100 PeV. Top: a simulated proton shower. Bot- 
tom: experimentally observed shower. 



distribution functions of the VEM signal having the same 
form as those used by the Haverah Park and Auger Obser- 
vatory. We have not studied the lateral distribution of sig- 
nal at distances smaller than 50 m because the lateral dis- 
tance and the uncertainty in the reconstructed core position 
become comparable to each other; this results in a large 
uncertainy in the signal at shorter distances. The LDF of 
both iron and proton showers becomes flatter with increas- 
ing zenith angle. 

We also see, as expected, the flatter lateral distribution of 
the simulated iron showers compared to that of the proton 
showers. The highest ratio of the iron to proton LDF occurs 
at distances more than 500 m from the shower core. We 
also show the distance at which the relative signal strengths 
of the signals from proton and iron showers are equal, the 
crossover radius depends significantly on the shower zenith 
angle and also on energy. For 100 PeV showers, it varies 
between 50-150 m for zenith angles as large as 25 degrees, 
and it is above 300 m at 45 degrees. For 10 PeV showers, 
not shown here, it varies between 130-180 m for zenith an- 
gles as large as 25 degrees, and it is above 400 m at 45 de- 
grees. For a mass independent energy reconstruction in this 
energy range, we will explore the possibihty of scahng the 
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1 Introduction 



The IceCube Neutrino Observatory started taking data in 
May 2011 with the complete array of 81 surface stations 
and 86 strings in the deep Antarctic ice [1]. The array grid 
is shown in Figure 1. In the period of June 2010-May 201 1, 
the surface air shower array, IceTop [2], was operated with 
73 stations (146 ice Cherenkov tanks) positioned on a trian- 
gular grid with a 125 m spacing. The IceCube detector had 
79 strings with 60 sensors on each string at depths between 
1450 m and 2450 m in the ice. We refer to this configura- 
tion as IceTop-73/IceCube-79. The detector collected data 
with 98% uptime during this period. 

IceCube, located at the geographic South Pole (altitude: 
2835 m), is at an optimum atmospheric depth of 680 g/cm^ 
where cosmic ray air showers in the PeV energy range are 
close to their shower maximum. In addition, with fast dig- 
ital electronics for signal processing and high resolution 
waveforms, IceTop is in a unique position to make a de- 
tailed measurement of the cosmic ray energy spectrum in 
this energy region in a few years. About 30% of showers 
trigger both detectors; these are called coincident events. 
The energy deposited along the kilometer long tracks of 
the penetrating muon bundles in IceCube, when combined 
with the energy deposited on the surface in IceTop, provide 
a mass composition sensitive measurement. 
In this paper we evaluate the performance of IceCube as a 
three dimensional Extensive Air Shower (EAS) array based 
on nine months of data from June 15, 2010 to March 15, 
2011. The data is split into an austral winter dataset (Jun 
15, 2010 - Nov 1, 2010) and an austral summer dataset 



IceTop tank 
^ DeepCore string 




Figure 1 : The surface map of IceCube in its completed config- 
uration. IceTop Stations 1, 7, 14, 22, 31, 79, 80 and 81 were not 
yet deployed in 2010. The circle with a 400 m radius shows the 
containment criterium for the reconstructed shower core position. 
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Figure 2; The differential rate of events in cosine zenith angle is 
shown in the top plot for IceTop-73 data and for coincident data. 
The dark (blue) and light (green) lines represent the austral winter 
and summer rates respectively. The ratio of winter to summer 
rates in the bottom plot shows the zenith angle dependence. 



(Nov 1, 2010 - Mar 15, 2011) to investigate the basic 
shower observables and their temporal behaviour due to at- 
mospheric changes. After minimal cuts the quality of the 
reconstructed events already reaches an accuracy to probe 
the inherent systematics left in data which may degrade the 
energy resolution and thus need to be studied further. 

This analysis uses the largest statistical data set which cov- 
ers two austral seasonal conditions. Earlier analyses [3, 4] 
did not account for atmospheric changes and could there- 
fore not combine data of more than one season. 
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Figure 3: The differential rate of log]^Q(Si25) for four zenith 
angle bands of equal solid angle in the region where IceTop is 
fully efficient. The plot quantifies the attenuation of showers with 
zenith angle. This attenuation is corrected when converting S125 
to primary energy. 
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Figure 4: The differential rate of logjQ(5i25) for the coincident 
IceTop-73/IceCube-79 events with a line fitted between and 0.5 
to guide the eye. 



2 Reconstruction of basic observables 

IceTop measures the Cherenkov light emitted by charged 
particles passing through the tanks. The timing informa- 
tion is used to reconstruct the aiTival direction, while the 
signal strength is used for the core position. The lateral 
distribution of the energy deposition by each shower is fit- 
ted to a function of the signal strength in Vertical Equiv- 
alent Muons (VEM) versus distance to shower core. The 
signal strength evaluated at 125 m from the reconstructed 
shower core, 5i25, is found by simulation studies to be 
sufficiently mass independent to primary energy for nearly 
vertical showers [5]. For coincident events also the Ice- 
Cube signals were used, together with a fixed reconstructed 
core position at the surface, to fit the direction using a func- 
tion which describes the muon bundle range-out [4]. This 
improves the angular resolution [6] and will be used to re- 
move unrelated IceTop/IceCube coincident events. 



2.1 Event selection 

Air shower events which trigger at least five IceTop sta- 
tions were reconstructed with the standard IceTop recon- 
struction procedure [7]. The events were selected if the 
reconstructed observables converged successfully, and the 
core location was reconstructed within a circle of 400 m 
radius from the IceTop array center as shown in Figure 1. 
On average, the reconstructed event rate was 2.55 Hz with 
a 12% vaiiability due to the barometric pressure changes. 

Events which tiigger both the IceTop and the IceCube ar- 
ray, but do not belong to the same air shower constitute 
an important background for the coincident data sample. 
These so-called random coincident events were cut based 
on the time difference between the signals in IceTop and 
the signals in IceCube. The zenith angle difference be- 
tween IceTop and IceCube direction reconstructions was 
also used to remove random coincident events. For a good 
energy loss reconstruction in the deep ice, the muon tracks 
which were not well contained by the detector volume (cor- 
ner clippers) were not used in the final sample. This basic 
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Figure 5: The weighted spectrum of log]^Q(Si25) for the IceTop- 
73 events with cos6 > 0.8 is shown for winter and summer 
months by (blue) squares and (green) stars respectively. An expo- 
nential function is fitted to the slope of the winter curve between 
0.7 and 1.1 and extrapolated to log]^o('S'i25) ~ 2.5 to guide the 
eye. 



set of cuts results in a high quality, well reconstructed data 
sample with an angular resolution below 1° and a core res- 
olution of about 10 m. The average IceTop-73/IceCube-79 
coincident event rate after this selection was 0.72 Hz with 
a 15% variability. 

In total 34.8 million events were left in the IceTop-73 data 
sample, while 9.5 million events remained in the coincident 
data sample (for a hvetime of 236.44 days). 

2.2 Basic shower observables and their temporal 
behaviour 

Figure 2 shows the differential rate in cosine zenith angle 
for both IceTop-73 events and IceTop-73/IceCube-79 coin- 
cident events. The requirement that the muon bundle com- 
ponent of the air shower must pass through the IceCube 
detector leads to a steeper distribution for the coincident 
events. The light (green) line is the rate in the austral sum- 
mer season and the dark (blue) line shows the rate in the 
austral winter season. The change in the atmospheric den- 
sity profile [8] from winter to summer as well as the snow 
accumulation of 21 cm on average during 2010 contribute 
to the observed difference in rates between the two seasons. 
The ratio quantifies the angular dependence of both atmo- 
spheric and snow effects. 

The differential rate of log]^Q(S'i25) is plotted for four 
zenith angle bands of equal solid angle in Figure 3 in the 
range where IceTop is fully efficient. The rate difference 
observed between the different zenith bands is due to the 
atmospheric attenuation. We could apply the classic Con- 
stant Intensity Cut analysis by calculating the attenuation 
for different zenith angles and evaluating the spectra at 
one zenith angle. In this analysis the IceTop-73 events 
were restricted to cos 6 > 0.8 to study the atmospheric 
effects within this narrow angular range. However, with 




Mean 0.44O4± 0.0045 



3 

•S, 0.09 
I 

0.085 
^ 0.08 
0.075 
0.07 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
log^^(S125/VEM) 



Figure 6: The knee region of the 5*125 spectrum. A shift of 0.04 
in log jg (5*125 ) is observed between the winter (blue circles) and 
summer (red triangles) spectra. 



high statistics it will be possible to use each zenith band as 
a separate dataset to reconstruct and analyze with different 
methods and further explore its mass composition sensitiv- 
ity. 

In Figure 4 the differential rates of Iog2o(<5'i25) for the co- 
incident IceTop-73/IceCube-79 events are presented. The 
break seen in slope of the spectra around log]^o('5'i25) = 0.4 
is the cosmic ray knee. 

To probe the details in the spectrum the rate weighted by 
for cos 9 > 0.8 is plotted in Figure 5. It is inter- 
esting to note the opposite behaviour of the high and low 
energy events in different atmospheric conditions, causing 
a change in the slope. The spectra show a hardening trend 
between 1.2 < log]^Q(S'i25) < 2.1. This systematically 
significant feature could not be traced back to any anomaly 
in data and is also seen by previous analyses [3, 4]. 

The knee position changes by 0.04 in logj^Q (5*125), from 
0.44 in summer to 0.40 in winter as shown in Figure 6. As 
the knee position should be at the same primary energy, 
the conversion from 5i25 to primary energy will need to 
take into account that the atmospheric variation and snow 
accumulation changes the detector response for the same 
primary energies. 

From an energy estimator, which was derived by simula- 
tions of proton showers, logio('5*i25) = 1 corresponds to 
about 10 PeV, logio( 6*125) = 2 to about 100 PeV, and 
logio(>5*i25) — 2.7 to about 1 EeV. From Figures 3 and 4, 
the measured rates indicate that about 150 events per year 
are measured by the three dimensional IceCube air shower 
array above 300 PeV and about 15 events above 1 EeV. 
There is also no indication of array saturation up to 1 EeV. 
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Fi gure 7; The average muon bundle energy loss measured 
by IceCube, shows a strong correlation with the shower size, re- 
constructed by IceTop, up to the highest cosmic ray energies. 

3 Sensitivity to primary composition 

When EAS propagate through the South Pole ice layer, 
only narrow 20 m) bundles of highly energetic muons 
(>500 GeV) survive to the depth of IceCube. They will 
lose energy mainly stochastically which then generates 
Cherenkov light detected by IceCube. Here the average en- 
ergy loss of muon bundles at the center of the IceCube de- 
tector volume is reconstructed using the measured charge 
signals and taking into account light propagation and the 
ice properties [9]. The energy loss of a muon bundle is a 
convolution of the muon multiplicity, the muon energy dis- 
tribution within the bundle and the single muon energy loss. 
Both the multiplicity and the energy distribution depend 
highly on the primary mass and energy. Thus, the muon 
bundle energy loss provides a composition sensitive ob- 
servable and can be determined with a resolution of about 
0.3 in logioof ^ [9]. 

In Figure 7, nine months of high quality coincident data 
is shown. There is a clear correlation between the IceTop 
energy sensitive observable, 5*125, and the energy loss up 
to the highest energies. Also, the shape of the distribution 
becomes narrower for higher 5*125 (larger showers). From 
simulations we know that proton and iron distributions be- 
come well separated in energy loss as function of primary 
energy [6]. The narrower distribution could therefore be re- 
lated to a change in composition and will be studied further 
by simulations. 

Events with a relatively low energy deposition in IceCube 
and no correlation with the IceTop shower size form the 
main background at this basic level of cuts as these are 
mainly caused by remaining random coincident events. 
However with minor cuts this is already reduced to the 
^1% level and will be further reduced in the next stage 
of the analysis. 



4 Summary/Outlook 

In this paper the performance of IceCube as a three 
dimensional cosmic ray detector was investigated using 
nine months of high quality data from June 15, 2010 to 
March 15, 2011 when the IceCube detector was in its 73- 
station/79-string configuration. We studied the behaviour 
of the basic observables in two distinct atmospheric con- 
ditions with the data split as austral winter and summer 
sets. We observed the change induced on the lateral de- 
velopment of showers by the atmospheric effects and snow 
accumulation, and a shift of 0.04 in logig (5*125) around the 
cosmic ray knee region between the two datasets. This ob- 
servable effect on the shower development will be studied 
in detail. We observed that the high and low energy show- 
ers are affected by the seasonal effects in opposite way, 
having a higher rate of high energy showers in the aus- 
tral summer atmosphere than in winter. The 5i25 spec- 
tra shows a systematically significant hardening between 
1.2 < log]^Q(S'i25) < 2.1. The cause of these effects, 
whether due to a change in shower maximum, change in 
the energy spectrum or a change in mass composition, will 
be the subject of further studies. 

The observed zenith angle dependence of EAS and the sea- 
sonal expansion of the atmosphere can be used for compo- 
sition studies. For coincident IceTop/IceCube air showers 
the composition sensitive ratio of muon bundle energy loss 
to shower size at the surface as well as the width of the 
energy loss distribution is being exploited. 

A few years of cosmic ray data will already provide enough 
statistics to accurately measure cosmic ray energy spectrum 
and composition up to 1 EeV by using multiple composi- 
tion and energy dependent shower observables through dif- 
ferent techniques such as neural networks [4]. 
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Abstract: Seven years of construction on the IceTop air shower array have culminated in the final detector setup 
of 162 ice-filled tanks. The tanks are paired as stations over an area of one square kilometer. A continuous and 
automatic procedure calibrates each tank via the extraction of the vertical equivalent muon peak in the tank charge 
spectrum. Over the years snow has drifted unevenly on top of the tanks. The overburden of snow influences the 
charge spectrum as the electromagnetic part of an air shower is attenuated more than the muonic part. The impact on 
individual tanks affects trigger rates and the air shower event structure. We will present tank response studies with Monte 
Carlo simulations and compare them to measured IceTop charge spectra in light of the varying thickness of the snow layer. 
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1 Introduction 

IceTop is the cosmic ray surface detector of the IceCube 
Neutrino Observatory [1], which is located at the South 
Pole and was completed in March 2011. The final setup 
consists of 162 tanks paired in 81 stations on a hexagonal 
grid with horizontal separations of about 125 m. The ar- 
ray covers an area of one square kilometer at an altitude 
of 2835 m above sea level, which is equivalent to an aver- 
age atmospheric depth of 680 g/cm^. The location together 
with the geometry and the applied trigger settings make 
the array sensitive to cosmic rays with energies between 
lOOTeVandlEeV. 

An IceTop tank acts as a calorimeter, converting energy 
deposited by charged relativistic particles into Cherenkov 
Ught. The reconstruction of air shower observables such as 
space angle (zenith and azimuth), core location (x and y) 
and shower size relies on a good understanding of the light 
signals in the individual IceTop tanks. A natural way to 
calibrate and simulate tanks with their individual responses 
is to use the reference signal - vertical equivalent muon 
(VEM) - that high energy vertical muons imprint in the 
charge spectra of each tank. 

A schematic of an IceTop tank is shown in Figure 1 . The 
tank is a cylindrical polyethylene vessel 186 cm in diam- 
eter and 130 cm in height covered with insulation. The 
walls are 6 mm thick and the inner part is covered with 
4 mm thick diffusely reflective liner of Zirconium' fused 
polyethylene. Each tank is initially filled up to a height 
of about 90 cm with purified water. Afterwards the water 



freezes via a controlled top-to-bottom procedure that lasts 
for about 40 to 50 days. To minimize impurities or contam- 
ination and cracks in the ice during the freezing process, 
a degasser system in combination with a circulation pump 
extracts air, swirls the water slowly around and removes ex- 
cess water. At the top of the ice volume, two Digital Opti- 
cal Modules (DOMs) are mounted to record the Cherenkov 
light pattern created by through-going and stopping parti- 
cles. A DOM [2] is an ensemble of a photomultiplier tube 
(PMT) and digitising, communication and operation elec- 
tronics inside a glass pressure sphere. In data taking mode 
the two PMTs inside one tank operate at a different voltage 
setting, enlarging the dynamic range in linear pulse charge 
assignment. Moreover, in case one DOM in a tank fails, the 
tank is not lost in further physics analysis. The ice surface 
and the two DOMs are covered with perlite, a volcanic ash, 
that acts as an insulator and prevents light from entering the 
ice volume and triggering the PMTs. 

The tanks are embedded in snow until the snow level equals 
the surrounding surface. So at installation there is only 
snow on the sides, but winds blow snow on top of the tanks. 
Each year some tanks accumulate up to 30 cm of snow on 
top while others are still in the original setup without a 
snow layer. This leads to an asymmetry in the behavior 
of the array and needs to be understood via data and Monte 
Carlo comparisons. 

1. Note that this is true for 150 tanks. The first 8 installed tanks 
in the year 2005 and 4 out of 16 tanks installed in the year 2010 
are different and contain a Tyvek bag that acts as a reflector. 
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Figure 1 : Schematic of an IceTop tank with its dimensions 
and its most important components. The trench where the 
tank is installed is backfilled with snow. Winds bring in 
snow that can accumulate on top of the tank. 

2 Calibration 

A good handle to study tank dependencies is to use the 
tank charge spectrum that is naturally generated by the very 
abundant low energy primary proton showers. Secondary 
muons within the 1 GeV to lOGeV energy region leave a 
very clear and distinguishable signal behind in the tanks. 
The technique was used to calibrate the water tanks in the 
Auger experiment [3], [4] and it is also in use in IceTop [5] 
since 2005. An example of a charge spectrum (expressed 
in number of photo-electrons (PE)) of one single tank can 
be seen in Figure 2. The data was recorded between the 
15th and the 22nd of January 2011. The overall spectrum 
is fit by the sum (solid) of a signal part (dashes) attributed 
to muons and a background part (dots) generated by elec- 
trons, positrons and gammas. The exact formula [6] reads: 



-2500 
£1 

E 



muon peak position = 130.28 PE - 
1 VEM position = 123.77 PE 
valley position = 91.01 PE 
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Figure 2: This is an experimentally measured charge spec- 
trum of DOM 61-61 running in calibration mode. The total 
fit (solid) is the sum of the muon (dashed) and electromag- 
netic (dotted) contributions. The three extra lines visualize 
the positions given by the legend. 
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Figure 3: This histogram represents the VEM values of all 
the 146 (high gain) DOMs running between the 15th and 
the 22nd of January 2011. The lines visualize the mean 
and the root mean square region. 
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where the first term in the summation describes the muon 
part: on its own a sum of a step function that describes cor- 
ner clipping muons and L(x,pi,p2) a Landau function with 
pi the most probable value and p2 a scale parameter that 
describes the bulk of through-going muons. So the signal 
part describes muons of all energies and incoming direc- 
tions. Muon telescope measurements have pointed out that 
vertical muons build up charges around 95 % of the overall 
muon peak. The last term in the equation describes the ex- 
ponential electromagnetic contribution. The three dashed 
vertical lines visualize the position of the valley, the muon 
peak and the 1 VEM unit. 

Since June 2009 all the tanks of the IceTop detector are cal- 
ibrated via an automated online procedure based on the fit 
described above. The muon calibration data is collected si- 
multaneously with the physics data and hence both DOMs 
are evaluated at their running^ gains. The calibration data is 
taken via a dedicated trigger that records every 8192nd hit 
per (high gain) DOM that is not in local coincidence with 
the other tank in the same station, i.e. there was no activity 
in the other tank within a short time window of 1 /is. 

Each IceTop tank responds in an individual way to air 
shower particles due to small differences in size, ice qual- 
ity, liner reflectivity, snow coverage and interfaces between 
the ice and DOM, the ice and the liners, the ice and the 
perlite and the ice and the degasser. The size of the ice 
block is important as more/fewer photons will be generated 
for longer/shorter tracks. Reflectivity and transmission will 



2. In the old method calibration was performed on a different 
set of data where both DOMs were adjusted to the same gain. 
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depend on the amount of impurities and cracks in the ice, 
the quality of the Hner and the interface between the ice 
and the surrounding material. The snow on top insulates 
the tank and so the PMTs in different tanks operate at dif- 
ferent temperatures and hence gains. Snow also blocks the 
low energy electrons, positrons and gammas and so the en- 
ergy threshold differs. Figure 3 shows the variety in PE per 
VEM for the 146 (high gain) DOMs operational between 
January 15th, 2011 and January 22nd, 2011. The average 
number is about 123 PE per VEM with a root mean square 
of 23 PE per VEM. The outliers with higher values belong 
to tanks Uned with the higher reflectivity Ty vek. 

3 Simulated VEM Spectrum 

The VEM spectrum is the result of the convolution of two 
components: the ground particles created in air showers 
and the detector response for each of the particles. The air 
showers in this study are produced via CORSIKA (v6900) 
[7]. The high energy hadronic interactions are treated via 
SIBYLL (2.1) [8] while the low energy interactions are 
taken care of by FLUKA (2008) [9]. The transition be- 
tween the models is set at center of mass energy of 80 GeV. 
The electromagnetic interactions are treated via EGS4 [10] 
parametrisations. The detector response to each particle 
is simulated via the official IceTray software [11] of the 
IceCube Collaboration. The IceTop tanks are modeled in 
a GEANT4 [12] environment that also tracks the ground 
level particles. The amount of Cherenkov photons that 
are generated inside the optical volume is used to assign 
the measured number of photo-electrons. The DOMs allo- 
cate a charge via an experimentally measured single photo- 
electron charge distribution. 

The air showers are generated according the following 
conditions: only protons are assumed; the zenith angle 
ranges from to 89 degrees and the azimuth angle from 
to 360 degrees. The azimuth angle is sampled uniformly 
while the zenith angle is sampled from a cos(0) sin(0) re- 
lation, the flat detector approach. The lowest triggering 
primary energy is given by the minimum amount of en- 
ergy needed to create a muon in the atmosphere that can 
possibly arrive at the surface. At the South Pole this en- 
ergy is about 3 GeV. The high energy limit is restricted 
by the caUbration trigger condition that deletes hits from 
the calibration stream when there is activity in the nearby 
tank. Simulations point out that the upper limit is situated 
around 30TeV. Moreover, Figure 4 shows that it is suffi- 
cient to simulate up to lOTeV as the integrals of the muon 
spectra saturate. The x-axis is expressed here in units of 
VEM/0.95, so that the muon peak is situated at one. In 
the threshold region [0;0.4] the highest energies contribute 
a bit more than in the muon peak region [0.8; 1.2]. This is 
also visible in Figure 5 where the integrals of the former 
spectra are shown as a function of the maximum energy in 
the dataset. Between lOTeV and 30TeV the integral in- 
creases about one to two percent, see therefore the solid 
line. The other markers give the integrals for the intervals 



10TeV 

2TeV 




Charge (VEM/0.95) 

Figure 4: Different simulated muon spectra (weighted to 
E^^-^) where the maximum energy that has been simulated 
changes. From bottom to top the maximum energy goes up 
between 100 GeV and 30TeV as written in the legend. 
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Figure 5: The same as Figure 4, but here the integrals of the 
muon spectra are plotted as a function of the maximum en- 
ergy Eniax used to build up the muon spectra. The solid line 
represents the whole integral. The dashed, dotted, dash- 
dotted lines are the integrals in following intervals [0;0.4], 
[0.4;0.8] and [0.8; 1.2]. 

[0;0.4] (squares), [0.4;0.8] (triangles) and [0.8;1.2] (inverse 
triangles). 

The VEM spectrum technique is used to calibrate the air 
shower simulations for high energy analysis. The simu- 
lated charge signals are calibrated in VEM units such that 
the simulated calibration muon peak is at the same posi- 
tion as in the corresponding experimental calibration data. 
The simulations perform well. All systematic checks com- 
pleted and listed in Table 1 exhibit a stability in the 1 VEM 
position within two to three percent. The first (second) col- 
umn gives the shifts of the 1 VEM position for simulations 
weighted to E^^ (E^^^). The checks included tests on 
stability of random numbers, hadronic interaction uncer- 
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Check 


E "'■^ 


E "'■^ 


different seeds 


0.2 % 


0.5 % 


July vs January atmosphere 


1.6% 


2.7 % 


SYBILL vs QGSJETOlc 


0.7 % 


2.9 % 


10 cm to 200 cm of snow 


0.4% 


0.7% 


different tanks/DOMs 


1.5% 





Table 1: Results of systematic checks for the shift of the 
muon peak and according 1 VEM position due to a vari- 
ety of effects listed. Two primary energy weightings are 
presented. 



tainty, daily variations of the atmosphere, the accumulation 
of snow and the differences of all (high gain) DOMs. 

4 Snow Dependency 

Over the years, snow drifts unevenly on top of the tanks. 
As can be seen from Figure 6, the overburden of snow in- 
fluences the charge spectrum as the electromagnetic part 
[0;0.75] of an air shower is attenuated more than the muonic 
part [0.75;1.25]. The ratio (S/B) between integrals of the 
muonic (M(h)) and the electromagnetic fit (EM(h)) in Fig- 
ure 2 and expressed in the following formula: 

^2.0VEM ^2.0VEM 

{S/B)ih) = / M{h)/ / EM{h) , 

J0.3VEM J0.3VEM 

is a function of snow depth h. This is shown in Figure 7. 
The markers represent simulations where the primary en- 
ergy is weighted via different differential spectral indices. 
The dotted and dashed lines are extrapolations to indices of 
-2.6 and -2.7 respectively. Direct weighting to these more 
realistic values of the differential spectral indices was not 
possible due to low statitics and large fluctuations. In com- 
parison the experimental relation is given by the solid line 
which shows a bit stronger dependency on the thickness of 
the snow layer. 

5 Summary 

The VEM spectrum can be efficiently generated with an 
uncertainty of 2-3%. This leads to an absolute charge cal- 
ibration in simulation. The snow dependency of the signal 
to background (muon versus electromagnetic) ratio in ex- 
perimental data is slightly stronger than the one found in 
simulated data and has to be further investigated. 
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Figure 6: The four histograms represent the (E~^'^) simu- 
lations of DOM 61-61. The only difference is the amount 
of snow on top of the tank lids. From top to bottom 
the snow thickness amounts to 10 cm, 40 cm, 70 cm and 
100 cm. 
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Abstract: While the primary purpose of IceCube is the search for high-energy astrophysical neutrinos, the overwhelming 
majority of events is caused by downgoing cosmic-ray induced muons. This provides a high-statistics data set which can 
be used for both detector calibration and supplemental physics investigations. In this work, we present a method to 
identify TeV-scale catastrophic energy losses along muon tracks and its application to the separation of single high- 
energy muons from large-multiplicity bundles which dominate the event sample above the horizon at high energies. The 
information can be used to derive the single-muon energy spectrum at all zenith angles up to energies of hundreds of 
TeV. We demonstrate that our measurement is sensitive to a cutoff of the proton spectrum at the cosmic ray knee and 
potentially to the prompt lepton flux caused predominantly by decay of charmed hadrons in atmospheric CR interactions. 

Corresponding authors: Patrick Berghaus'^ (berghaus@icecube.wisc.edu), Chen Xu^ (chen@udel.edu) 

^Bartol Research Institute, University of Delaware, USA 

Keywords: Muon Spectrum, Stochastic Losses, Knee 



1 Introduction 

In recent work on high-energy atmospheric muon fluxes, it 
was pointed out that second-generation astrophysical neu- 
trino detectors such as IceCube should be able to substan- 
tially extend the energy range of muon energy spectrum 
measurements and hence address various important open 
problems in cosmic -ray physics. These include verifica- 
tion of hadronic interaction models, primary composition 
at Pe V energies and prompt contributions to the lepton flux 
[1]. 

From an experimental perspective, the main challenge in 
the measurement of the muon energy spectrum is limited 
spatial resolution due to the sparse instrumentation of the 
detector. With a minimum separation of 17 meters between 
Digital Optical Modules (DOMs) in IceCube, resolution of 
individual muons within a shower is usually impossible. 
Previous measurements [2, 3] took advantage of the shorter 
range of large-multiplicity muon bundles compared to in- 
dividual high-energy muons during passage through matter 
to measure the muon spectrum at large slant depths near 
the horizon. The limited angular resolution of the detec- 
tor and background from up-going atmospheric neutrinos 
mean that the range of this technique is limited to values 
up to approximately 100 TeV. Furthermore, the muon spec- 
trum can in this way only be measured in an angular region 
where the relative contribution from prompt interactions is 
lowest, excluding investigation of this flux component for 
practical purposes. 
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Figure 1 : Probability for a stochastic loss above a given en- 
ergy to occur within the detector volume, for three typical 
muon energy ranges. The simulation was performed using 
MMC [5]. 

In order to take advantage of the full potential of a large- 
volume detector, it is therefore necessary to develop a 
method that allows separation of showers with exception- 
ally highly energetic muons from the dominant background 
of bundles containing low-energy muons of higher mul- 
tiplicity. It has long been proposed to make use of the 
fact that the stochasticity of energy deposition along muon 
tracks increases with energy [4]. In this analysis, by iden- 
tifying muon track segments with unusually strong photon 
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■ Dust Layer (-lOOni to -50m) 
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■ Ordinary Ice (200m to 250m) 
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Figure 2: Left: Average number of photo-electrons registered in IceCube DOMs as a function of perpendicular distance 
from reconstructed muon tracks for three different depth bins. Values are averaged over events with a total charge Qtot > 
1000 p.e. Right: Value of scattering parameter X^xp from fitting distributions on the left to Eq. 1 in dependence of vertical 
depth. The value of zoom — corresponds to the center of the IceCube array, located 1950 meters below the surface 
of the ice. The effective scattering length in two different ice models [6] is included for comparison. It is important to 
note that the physical meaning of X^ff and X^xp is not identical. Fine structures in the depth dependence are inherently 
smeared out in the experimentally derived parameter. 



emission, the amplitude of individual stochastic losses can 
be used to infer the most likely energy of the parent muon, 
as illustrated in Figure 1 . 

2 Energy Estimation Method 

The IceCube detector array records Cherenkov photons 
emitted by relativistic particles during passage through ice. 
The amount of charge in the optical modules provides 
calorimetric information that can be used to calculate the 
energy released in the event. 

Previous energy reconstructions relied on prior assump- 
tions about ice properties whose accuracy was necessarily 
limited. The method described here is based exclusively 
on experimental observations. Its is based on two simple 
assumptions: 

• Down-going tracks in IceCube consist mainly of 
muon bundles dominated by minimum-ionizing 
muons which lose their energy smoothly and contin- 
uously. 

• The total number of emitted Cherenkov photons is 
in good approximation proportional to the deposited 
energy. 

As illustrated in Figure 2, the total amount of light regis- 
tered in the DOMs can then be approximated by the empir- 
ical function 

dE{x = XQ) \ \ /d(xo)-25\ 

where ^^^j^ is the energy loss at vertical depth x, Q is the 
total charge in each DOM and d is the perpendicular dis- 
tance between a DOM and the reconstructed muon track. 



This data-derived model only relies on a single scattering 
parameter Xexp, which varies in dependence of the local 
ice properties. From the charge in each DOM, it is possi- 
ble to calculate a value that is approximately proportional 
to the energy deposited in its vicinity. Monte-Carlo simu- 
lated data are only required to determine the proportional- 
ity factor in Eq. 1 and for verification. Assuming point- 
like isotropic emitters on the track, clusters of exception- 
ally high values can be used to identify energetic stochastic 
energy losses. An example can be seen in Figure 3. 

There are two important additional benefits that arise natu- 
rally from this technique. The first is the reduction of biases 
resulting from uncertainties about exact ice properties, al- 
lowing for detailed investigation of simulation biases. The 
other is the availability of information about differential en- 
ergy deposition, permitting the construction of more accu- 
rate event energy estimators. 

3 Analysis Procedure 

This analysis was based on data from IceCube in its 59- 
string configuration (IC59). Only events with a total of 
more than one thousand photo-electrons registered in the 
detector were considered. This selection provides a sam- 
ple of events at all zenith angles which is not biased by the 
quality criteria used to select physics filter streams. 

The analysis procedure itself mainly consists of two rela- 
tively simple steps. First, a sample of muon tracks with 
strong stochastic losses inside the detector volume were 
identified. Then, the energy of the strongest single energy 
deposition was reconstructed and its spectrum compared to 
various benchmark models. 

The principal conditions that were imposed in the selection 
are: 
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Figure 3: Example for massive stochastic loss event found in IC59 data. The energy deposition estimate shown on the 
right hand side is consistent with an 80 TeV cascade. The x-axis corresponds to the distance along the track, measured 
from the surface of the ice. 
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Figure 4: Left: Muon multiplicity in cosmic-ray showers at entry into detector volume before and after selection cuts, 
simulated using CORSIKA/SIBYLL. The analysis favors muon bundles of moderate multiplicities while reducing both 
high-multiplicity bundles and low-energy tracks with high recorded DOM charges due to unusual topologies. Right: 
Effect of selection cuts in dependence of true maximum stochastic loss energy. The effective threshold energy lies around 
10 TeV. 



• Presence of a well-reconstructed muon track with a 
length of at least 600 meters inside the detector vol- 
ume. 

• A peak energy loss value in a 50-meter bin along the 
track which exceeds the median differential energy 
loss in the event by a factor of 10. 

• Location of peak energy loss within 150 meters from 
at least 50 DOMs, assuring containment within the 
main detector volume. 

Figure 4 shows the cumulative effect of the cuts. 

So far, the effective threshold energy for stochastic losses 
lies around 10 TeV. It is in principle possible to reduce this 
threshold further by relaxing the requirement on the total 
recorded DOM charge. 



4 Result 

To estimate the sensitivity of the method, the measured en- 
ergy spectrum of stochastic losses was compared to simu- 
lations based on various primary composition models [7]: 

• Constant Composition: The spectrum of each pri- 
mary nucleus type changes by the same amount at an 
energy corresponding to the cosmic-ray knee. 

• Rigidity-Dependent Cutoff: Each component of 
the cosmic ray flux is subject to a cutoff at an energy 
that is proportional to the charge of the nucleus. 

• Mass-Dependent Cutoff: Similar to the previous 
model, but assuming a cutoff dependent on the nu- 
clear mass. The main distinction here is a sharper 
transition, resulting in a stronger cutoff signature. 

• No Knee: This additional model was included for 
purposes of illustration only, assuming unbroken 
spectra for each component and therefore completely 
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eliminating the knee. It should be noted that this as- 
sumption is strongly disfavored by previous experi- 
mental results. 

Figure 5 shows the preliminary result, based on a subsam- 
ple of 10 percent of the data taken with IC59. For models 
incorporating a cutoff in the individual primary spectra, a 
corresponding effect can clearly be discerned in the simu- 
lated muon-induced stochastic losses. The limited exper- 
imental statistics do not yet allow any definite statement 
about preferred models, even while neglecting systematic 
uncertainties. 




10' 

























































1 










m 




Ww 








^mWW 




• 


• 




1 


PRELIMINARY 


T 


: -1- 
-•- 










— • — Data 

— ■ — Constant Composition 

— * — Mass-Dependent Ay 

— ¥ — Rigidity-Dependent A Y 

— e— No Knee ■ ■ 

1 1 1 1 ■VI- 





Ecascfi, [GeV] 



Figure 5: Comparison of reconstructed maximum stochas- 
tic loss energy for simulated and experimental data 
summed over all zenith angles. All simulated curves were 
obtained by reweighting the same data set generated with 
CORSIKA. Experimental data correspond to a subset of ten 
percent of IC59 events. The error band is purely statistical. 



5 Conclusion and Outlook 

As can be seen from comparing Figures 5 and 6, the highest 
muon energy found in the analyzed sample is most likely 
located around 500 TeV, with slight variations depending 
on the exact primary spectrum. This means that despite the 
fact that the spectrum falls off almost as the fourth power 
of the particle energy, the range of the measurement can 
be extended by roughly an order of magnitude compared 
to earlier experiments. The method presented here further- 
more allows the measurement of cosmic-ray induced muon 
spectra independently of the zenith angle, introducing an 
entirely new degree of freedom. 

Consequently, it will be possible to address a variety of new 
physics issues that were previously inaccessible. It has al- 
ready been demonstrated that the behavior of cosmic rays 
around the knee can now be probed directly. A separate but 
equally important question is the contribution of prompt 
production processes to lepton fluxes at high energies. In 
optimistic models, the charm-induced component becomes 
dominant around 300 TeV at near-vertical angles [8] and 
would therefore be directly measurable. If on the other 



Figure 6: Relation between reconstructed cascade and 
muon surface energy obtained from Monte-Carlo simula- 
tion, assuming a poly-gonato primary spectrum and neglec- 
tic prompt muons. The change in slope around 10 TeV can 
be attributed to threshold effects. 

hand the prompt muon flux were low enough to be dom- 
inated by decays of unflavored light mesons [9], it should 
be possible to set a stringent limit that would strongly con- 
strain the atmospheric neutrino background in searches for 
astrophysical sources. 

It should also be emphasized that the method to distinguish 
unusually energetic muons from large-multiplicity back- 
ground events is readily adaptable for analyses targeting 
neutrino fluxes at energies in the PeV range. 
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Abstract: Muons with a high transverse momentum (px) are produced in cosmic ray air showers mostly via semileptonic 
decay of heavy quarks and the decay of high pr kaons and pions in jets. These high pt muons have a large lateral 
separation from the shower core and accompanying muon bundle. IceCube, a kilometer-scale neutrino telescope 
consisting of an array of photodetectors buried in the ice of the South Pole and a surface air shower array, is well suited 
for the detection of high pr muons. The surface shower array can determine the energy, core location and direction of 
the cosmic ray air shower while the in-ice array can do the same for the high pr muon. This makes it possible to measure 
the cosmic ray muon lateral separation distribution at distances greater than 200 meters. The preliminary results from 
analysis of data from 25% of the full IceCube detector will be presented. 
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1 Introduction 

IceCube, a kilometer-scale neutrino telescope, is well 
suited for measuring the lateral separation of muons in cos- 
mic ray air showers. Completed in December 2010, it con- 
sists of a 1 km^ array of optical sensors (digital optical 
modules, or DOMs) buried deep in the ice of the South Pole 
and a 1 km^ surface air shower array called IceTop. IceTop 
has an energy threshold of 300 TeV and can reconstruct 
the direction of showers with energies above 1 PeV within 
'^1.5° and locate the shower core with an accuracy of 
^10 m [1]. The in-ice DOMs (here referred to as IceCube) 
are buried in the ice 1450 m under IceTop on kilometer- 
long strings of 60 DOMs with an intra-DOM spacing of 
17 m. IceCube can reconstruct high quality tracks of high 
energy muons with < 1° accuracy. IceTop and IceCube can 
be used together to select cosmic ray events with a muon 
with a minimum lateral separation of ~200 m. The mea- 
surement of the lateral separation distribution of muons in 
air showers provides a valuable check on air shower sim- 
ulation models and can be used as an independent method 
for determining the cosmic ray composition [2]. 

The most common source of muons with a large lateral sep- 
aration is muons with a high transverse momentum (pt) 
primarialy from charm and bottom mesons and j ets of high 
Pt partons [2]. The transverse separation is given by: 



where is the energy of the high pT muon, and h is the 
interaction height of the shower, here taken as an average 
value of 25 km. The interaction height loosely depends on 
the composition and a full treatment of this is planned in the 
future. Using a separation of 200 m (75 meters more than 
the string separation) as a rough threshold for the two-track 
resolution distance of the high pT muon from the shower 
core gives a minimum resolvable pt of 8 GeV/c for a 1 TeV 
muon. Rough calculations predict on the order of tens of 
high Pt muon events in the studied data sample [2]. 
The combined acceptance for cosmic ray air showers that 
pass through both IceTop and IceCube is 0.3 km-^sr for the 
full 86-string IceCube array [3]. By the end of the austral 
summer of 2006/2007, 22 IceCube strings and 26 IceTop 
tanks had been deployed. The combined acceptance for 
showers that trigger both IceTop-26 and IceCube-22 is 0.09 
km^sr In 2007 the discriminator threshold settings for Ice- 
Top were changed partway through the year to sub-optimal 
values, so this analysis discards the data taken during this 
period leaving 114 days of livetime. A search has been 
conducted in this data for cosmic ray events with a muon 
with a large lateral separation. 

The underground muon detector MACRO has previously 
measured the separation between muons in air showers for 
shower energies ranging roughly from 10^ GeV to 10^ GeV 
[4]. MACRO measured muon pair separations out to a dis- 
tance of about 65 meters. They verified the hnear relation- 
ship between pt and separation shown in Eq. 1 (with a 
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Figure 1 : Zenith angle resolution of the high pT reconstruc- 
tion algorithm. The sigma are the results of Gaussian fits 
to the distributions. 



The resolution is worse for the high pT muon because 
fewer DOMs are hit. While high pT muons with a greater 
separation are much easier to resolve with the two track al- 
gorithm, they also tend to be lower energy (see Eq. 1) and 
hit fewer DOMs. The average number of DOMs hit by the 
high Pt muon is 15, compared to 90 for the muon bun- 
dle. Additionally, because high px muons have fewer hit 
DOMs, a DOM hit by the muon bundle that is incorrectly 
placed in the high pT group has a much larger effect on 
the reconstruction of the high pT track (it can also degrade 
the bundle resolution, but to a lesser extent). These factors 
degrade the resolution of the high pT track direction and 
timing. The spatial resolution of the high pT muon track 
(as measured by the difference between the reconstructed 
track and the true track at the point of closest approach to 
the detector) is 40 m in x, y, and z. 



small offset due to multiple scattering of the muons) out to 
momenta up to 1.2 GeV/c. 

2 Reconstruction Methods 

High Pt muons appear as a separate track coincident in 
time and parallel with the track from the central core of low 
Pt muons. Generally the bundle of low pT muons leaves 
more light in the detector than the high pT muon. 

Current reconstruction algorithms in IceCube are designed 
to reconstruct single tracks. In order to reconstruct high 
Pt double-track events the hit DOMs (i.e. DOMs that de- 
tect at least one photon) are attributed to either the muon 
bundle or high pT muon based on their timing and posi- 
tion. The muon bundle is reconstructed first by the IceTop 
surface array. This reconstruction serves as a seed for a 
likelihood-based reconstruction using the IceCube DOMs. 
Only IceCube DOMs within 90 m of the seed track are 
used for bundle reconstruction. The high pT muon hits 
are selected relative to the reconstructed bundle track. The 
high Pt muon arrives at the same time but laterally sepa- 
rated from the bundle, so its hits will have a negative time 
residual relative to the bundle reconstruction. Additionally, 
high Pt muon hits are required to be at least 90 m from 
the reconstructed muon bundle track. The high pT muon 
hit track is reconstructed using a likelihood reconstruction 
with a downgoing hypothesis. 

Figure 1 shows the performance of this procedure after re- 
construction quality selection criteria have been applied. 
The zenith angle resolution for groups determined using the 
true simulation information (black, solid lines) is compared 
to the resolution for groups determined using the splitting 
algorithm (red, dashed lines) for the muon bundle (top) and 
high pt muon (bottom). Roughly 50% of the events fail to 
reconstruct because there are not enough DOMs in one of 
the groups. 

These reconstruction algorithms achieve a zenith angle 
(space angle) resolution of 0.6° (1.3°) for the muon bundle 
and 2.6° (8.2°) for a high pT muon separated by 200 m. 



3 Signal and Background Separation 

Since high pT muons occur in only a fraction of simulated 
showers, a toy model based on CORSIKA [5] proton show- 
ers was used to model the signal. A single muon is inserted 
into an existing CORSIKA event containing a muon bun- 
dle from an air shower. This modified shower is then run 
through the standard IceCube propagation, detector sim- 
ulation, and reconstruction routines. Simulations insert a 
muon with energy of 1 TeV separated 100, 150, 200, and 
400 m from the shower core. 

Cosmic ray air showers that do not generate a high pt 
muon (called 'single muon bundles') are a background to 
this search. Since they generate only one track in the ar- 
ray, these events are mostly eliminated by requiring there 
be two well-reconstructed tracks in the IceCube detector. 

The IceCube 22-string configuration is large enough that 
the rate of simultaneous events from cosmic rays is signif- 
icant. Muon bundles from two (or more) air showers can 
strike the array within the 10 fis event window, producing 
two separated tracks. These so-called double-coincident 
events are another background for air showers with high 
Pt muons. Since these double-coincident events are un- 
corrected in direction and time, requiring that both recon- 
structed tracks be parallel (within 15° of each other) and 
occur within ±600 ns can eliminate most of these events. 
However, an irreducible background remains from double- 
coincident events that happen to come from roughly the 
same direction and time. The rate of double-coincident 
events can be estimated by looking at the off-time rate of 
events (i.e. events with tracks that occur more than 600 ns 
apart) in the data that pass all other selection criteria. 

A number of selection criteria are applied to separate events 
with high Pt muons. The events are required to trigger 
at least 6 DOMs in IceTop and at least 8 DOMs in Ice- 
Cube. The events are also required to have high quality 
two-track reconstructions. This includes requiring that the 
tracks have at least one hit DOM with a time residual of 
less than 15 ns (a 'direct hit') and that tracks pass within the 
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detector fiducial volume. Further reduction of single muon 
bundles is done based on the differences in event topologies 
and timing. For instance, single muon bundles are well- 
reconstructed by a single track hypothesis, while the high 
Pt muon events are not. Figure 2 shows the negative log 
of the reduced likelihood of a single track reconstruction 
for single muon bundles, and showers with an inserted 8 
and 16 GeV/c pT, 1 TeV muon (separation of 200 m and 
400 m from the shower core, assuming an average inter- 
action height of 25 km). Well-reconstructed events have a 
lower likelihood value on this plot. For large separations, 
this variable separates single muon bundles from showers 
which contain a high pT muon. This analysis retains events 
with a likelihood greater than 7.5. Next we require that the 
high Pt muon track be a robust track in the detector with at 
least 6 direct hits. Finally, the remaining background single 
muon bundle events are removed by requiring the perpen- 
dicular separation between tracks to be at least 160 m. The 
number of events passing each type of selection criteria for 
data and simulated background and signal are shown in Ta- 
ble 1. 
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Figure 2: Negative log of the reduced likelihood for the 
single track reconstruction 



4 Results 

After applying all selection criteria 53 events remain in 1 14 
days of data. No events remain in the simulated back- 
ground. One event that passes all the selection criteria is 
shown in Fig. 3. The reconstructed high pT muon and bun- 
dle track have a perpendicular separation of 207 m at the 
center of the detector The two tracks arrive within 470 ns 
and 3.3° of each other 

4.1 Purity of Final Sample 

Simulation of single muon bundle background events is 
too computationally intensive to accumulate large statis- 
tics. Although no simulated single muon bundle event 
passed all the selection criteria, the possibility exists that 
some single muon bundle background events could survive 
in the final data sample. In order to estimate the purity of 



Figure 3: Candidate shower with a high px muon. The 
cosmic ray bundle is on the left and the high px muon is on 
the right. 



the final sample a study has been conducted of simulated 
events which pass the reconstruction quality and double- 
coincident seleciton criteria (without applying the criteria 
designed to remove single track events). At this cut level, 
the only events left in the single muon bundle background 
are events which are incorrectly split into two tracks by the 
splitting algorithm. Two previously unused variables were 
developed that focused on studying how well the bundle fit 
described the high px muon hits and vice versa. 

The first variable is the standard deviation of the time resid- 
ual relative to the reconstructed bundle track (abundie)- 
Only hits belonging to the high px muon reconstructed 
track (time residual < 100 ns and distance from the bun- 
dle > 90 m) were used to calculate cTi^n^idiQ. In the single 
muon bundle background abundie is small because the bun- 
dle reconstruction is equally good at describing both sets of 
hits. Conversely, events with high px muons have a later- 
ally separated track that is causally disconnected from the 
bundle track. This leads to larger values of abundie ■ 
The second variable is the mean of the time residuals rel- 
ative to the reconstructed high px muon track. Only hits 
within 90 m of the bundle track are used for this calcula- 
tion. This variable uses the fact that the parallel but lat- 
erally separated high px muon track have a negative time 
residual relative to the bundle track. Single muon bundles 
incorrectly split into two groups tends to have time residu- 
als that are closer to zero because both groups are causally 
consistent with the bundle hypothesis. 

These two variables are shown in Fig. 4. The back- 
ground single muon bundles and high px muon simulation 
are at a relaxed cut level with only reconstruction quality 
and double-coincidence selection criteria applied, while the 
data points are after all selection criteria have been applied. 
The dashed grey line is a fit to the mean of the background 
single muon bundle points. To estimate the purity of the 
data sample, this line has been moved upwards until every 
simulated single muon bundle event is below it. Thirty-one 
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Cut 


Data 


Simulation 


Simulated Signal (%) 


IceTop and IceCube Trigger 


1.35 X 10'' 


1.47 X 10 '' 


100% 


Bundle and High -pT Reconstructions Successful 


4.59 X 10^ 


3.04 X 10^ 


49% 


Double-Coincidence Cuts 


1.16x 10^ 


1.35 X 10^ 


25% 


Reconstruction Quality Cuts 


2.57 X 10^ 


1.64 X 10^ 


13% 


Single Track Cuts 


53 





8% 



Table 1: Number of events in 114 days for data and simulated background. The percentage passing rates for simulated 
high 'Pt muons with a lateral separation of 200 m (8 GeV/c) are shown as well. 
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Figure 4: Standard deviation versus mean of time residuals. 
Note that the selection criteria have been relaxed for the 
simulated events. 



of the fifty-three data events lie above the line, giving a pre- 
liminary estimation for the purity of the final event sample 
of at least 58%. 



These events have been scaled down to the rate of 3 events 
in the on-time window. 
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Figure 5: Preliminary perpendicular separation between 
bundle and high muon reconstructions at the center of 
the detector. 



4.2 Estimation of Double-Coincident Rates 

As mentioned in Section 3, an irreducible background re- 
mains from double-coincident cosmic rays. The number of 
these events in the on-time window (±600 ns) can be es- 
timated by counting the number of events with tracks that 
occur with a time difference greater than 600 ns after ap- 
plying all other selection criteria. In the final data sample, 
eight events were found in this off-time window. To be con- 
servative, two events with a time difference greater than 
20,000 ns are discarded. The remaining six events were 
spread over a time range of 2260 ns, giving an expectation 
that 3 of the 53 events in the ±600 ns on-time window are 
due to double-coincident cosmic rays. 

4.3 Lateral Distribution of Muons 

Figure 5 shows the preliminary perpendicular separation 
between the bundle and high muon reconstructions at 
the center of the detector The black solid lines are the 
53 data events that passed all selection criteria. Prelimi- 
nary estimations indicate that this sample is 58% pure. The 
dashed red lines are the data events that passed all selec- 
tion criteria, but had time differences greater than ±600 ns. 



5 Conclusions 

A search for cosmic ray events with laterally separated 
muons has been conducted in 114 days with 25% of the 
full IceCube detector. The fraction of single muon bun- 
dle events remaining in the final sample has been con- 
servatively estimated to be 42% and the number of back- 
ground double-coincident events has been shown to be very 
small. Further searches with IceCube data with larger in- 
strumented volume will have important implications for 
cosmic ray composition and air shower simulation. 
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Abstract: IceCube and its surface array IceTop can be used to search for PeV gamma rays. At this energy, pair production 
on the cosmic microwave background limits the gamma ray horizon to a few tens of kiloparsecs. Therefore, the gamma 
rays have to originate from within the Galaxy, either from the sites of very energetic accelerators or from the interstellar 
medium and dense molecular clouds, which act as a target for the local cosmic ray flux. 

IceCube can efficiently distinguish PeV gamma rays from the background of cosmic rays. For showers within a limited 
zenith angle range, the in-ice component of IceCube can detect the Cherenkov radiation released by muons that penetrate 
kilometers deep into the ice while the surface component, IceTop, can measure particles from the same air shower. 
Gamma-ray air showers have a much lower muon content than cosmic ray air showers of the same energy. IceTop can be 
used to search for gamma-rays by selecting those showers that lack a signal from a muon bundle in the deep ice. 

We investigate the possibility of detecting gamma-ray air showers with IceCube and present results of one year of data, 
taken in the 2008/2009 season when the detector consisted of 40 strings and 40 surface stations. We derive an upper limit 
on the ratio of gamma showers to cosmic rays showers of 8.1-10~* above 1.2 PeV in a region of 10 degrees around the 
visible part of the Galactic plane. The projected gamma-ray sensitivity of full IceCube, completed at the end of last year, 
is also studied. 
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1 Introduction 

Gaimna rays are an important tool for studying the cosmos; 
unlike cosmic-rays (CRs), they point back to their points 
of origin. Air Cherenkov telescopes have detected many 
high-energy {E > 1 TeV) gamma-ray sources. Within 
our galaxy, gamma-rays have been observed coming from 
several types of sources: supernova remnants, pulsar wind 
nebulae, binary systems, and the Galactic center. Extra- 
galactic sources include Active Galactic Nuclei (AGNs), 
and other objects containing supermassive black holes. 
Surface air shower arrays like Milagro have performed all- 
sky searches for TeV gamma-rays. Although these detec- 
tors are less sensitive to point sources than Air Cherenkov 
telescopes, they have identified several diffuse sources [1]. 

At higher energies, extra-galactic sources are unlikely to be 
visible, because more energetic photons are likely to inter- 
act with cosmic microwave background radiation (CMBR) 
photons, and infrared starlight from early galaxies, produc- 
ing e+e" pairs. At 1 PeV, for example, near the peak of the 
cross-section, photons are limited to a range of about 25 
kiloparsecs. It is unknown whether Galactic accelerators 
exist that can produce gamma rays of such high energy, but 



an expected flux results from the interaction of CRs with 
the interstellar medium (ISM). 

To date, the best statistics on photons with energies above 
about 300 TeV come from CASA-MIA, which has set a 
limit on the fraction of photons in the cosmic-ray flux of 
10~^, at energies above 600 TeV [2]. They also set a limit 
on the relative flux of photons coming from within 5 ° of 
the galactic disk of less than 2.4 x 10~^ [3]. This is close 
to the theoretical expectation due to CR interactions with 
the ISM. For a Northern hemisphere site like CASA-MIA, 
the average theoretical prediction is 2 x 10 ~^ [4] . At some- 
what lower energies, gamma-ray emission from the Galac- 
tic plane may already have been seen [5]. 

2 Detection Principle 

We select a sample of air showers that are successfully re- 
constructed by IceTop and have a shower axis that passes 
through the instrumented area of IceCube deep in the ice. 
Investigating the muon component of the air showers in the 
ice allows for CR composition studies. In this analysis, 
we create a sample of gamma-ray candidates by selecting 
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Figure 1 : Gamma-ray field of view for IceTop/IceCube in 
Galactic Coordinates (black circle). The rectangular area 
is the region for which CASA-MIA [3] has set upper hmits 
on the Galactic diffuse photon flux. The part of the Galac- 
tic Plane covered by IceCube is smaller but closer to the 
Galactic Center. 



only those air showers which do not produce any signal 
in the deep ice. The energy threshold of the air showers 
must be chosen sufficiently high to reduce the background 
of muon-poor CR showers. At 1 PeV, less than 0.1% of 
all proton showers do not contain muons with an energy 
over 500 GeV, which is approximately needed to reach Ice- 
Cube when traveling vertically. However, only 16% of all 
gamma-ray showers will produce a signal in IceCube at that 
energy. At lower energies, the CR background becomes too 
large to find a possible gamma-ray signal . Moreover, the 
shower size becomes to small to trigger a sufficient amount 
of IceTop tanks for high quality angular reconstruction. 

The requirement of IceTop-IceCube coincidence confines 
the maximum zenith angle to ^30 degrees. Since IceCube 
is located on the geographic South Pole, this means that 
our Field of View (FOV) is constrained roughly to the dec- 
lination range of -60 to -90 degrees, which includes part 
of the Galactic Plane (see Fig. 1). The Magellanic clouds 
also lie within the FOV, but at distance of ~50 kpc, a possi- 
ble gamma-ray flux at ^ 1 PeV from the Large Magellanic 
Cloud is reduced by a factor ~ 100. 

Within our candidate sample we search for a local increase 
in the flux from a point source or the Galactic Plane. We 
make use of the unique geographical position of IceCube: 
as the Earth rotates the declination remains a fixed func- 
tion of the zenith angle, while the Right Ascension (RA) 
rotates with respect to the azimuth angle. A diffuse flux 
of CR showers therefore produces a flat distribution in RA 
(whereas the azimuth distribution is a complex function of 
the detector shape). Therefore, a clustering of events in RA 
can be attributed to a gamma-ray component in the flux. 



3 IC40 Analysis 

Between April 1 2008 and June 1 2009, IceCube has taken 
data with a configuration of 40 strings and 40 surface sta- 
tions (IC40). There are several trigger conditions based 
on different signal topologies. For this analysis we use 
the 8 station surface trigger, which requires a signal above 
threshold in both tanks of at least 8 IceTop stations. Al- 
though an additional signal in IceCube is not required for 
this trigger, hits in the deep detector within a time window 
of 10 fxs before and after the surface trigger are stored. 

The air shower parameters are reconstructed with a series 
of likelihood maximization methods [7] . The arrival direc- 
tion is determined with a resolution of 1.5°. 
In the selection of photon shower candidates from the data 
sample we distinguish two steps: level three (L3) and level 
four (L4). At L3, two parameters are used to constrain the 
geometry and ensure the shower axis passes through the 
instrumented volume of IceCube . The IceTop containment 
parameter Cit is a measure of how centralized the core 
location is in IceTop. Cit = x means that the core would 
have been on the edge of the array if the array would be 
X times its actual size. The string distance parameter dgtr 
is the distance between the point where the shower reaches 
the depth of the first level of DOMs and the closest inner 
string. Inner string, in this sense, means a string which is 
not on the border of the IC40 configuration. IC40 contains 
17 inner strings. The complete set of L3 cuts is composed 
of these geometrical parameters, the reconstructed shower 
energy E^-^^, and reconstruction quaUty parameters. The 
L4 stage imposes only one extra criterion to the event: there 
should be no hits in IceCube. 

Simulations show that 16% of the gamma-ray showers is 
excluded by the L4 cut. The amount of CR showers still 
present at L4 is not simulated. Instead, this value is de- 
termined from data, where it is assumed that CRs domi- 
nate the sample in the off- source region (i.e. away from the 
Galactic Plane). 

3.1 Galactic Plane Test 

The IC40 data set consists of 368 days of combined Ice- 
Cube/IceTop measurements. The data from August is used 
as a burn sample, which means that it is used to tune the pa- 
rameters of the analysis, and then discarded. The remain- 
ing data is used for the actual analysis .The cut parameters 
of this analysis were optimized to maximize the chance for 
discovery, because the potential of this analysis to constrain 
the gamma ray flux would not be competitive with already 
existing limits. A discovery above existing limits is in prin- 
ciple possible, since our FOV covers a part of the Galactic 
Plane not seen by previous searches, close to the Galactic 
center. 

There are three cut parameters that are optimized by using 
the bum sample: E^j^^ , Cit and d^u- . This is done by scan- 
ning through all combinations of parameters within a cer- 
tain range: 600 TeV< E^ < 2 PeV in steps of 100 TeV, 
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Right Ascension (radians) 

Figure 2: The 268 candidate gamma-ray events superimposed on HI column densities based on [6] 



0.5 < CiT < 1.0 in steps of 0.1 and 50 m < d^tr < 90 in 
steps of 10 m. For each combination the number of events 
Ns is counted which are located in the source region, de- 
fined as within 10 degrees of the Galactic Plane. Then, 
the events are scrambled multiple times by randomizing the 
RA of each data point. The quality of the combination of 
cut parameters is given by the fraction of scrambled data 
sets that produce the same or higher amount of events in 
the source region as the original data set. The combina- 
tion i^jjjjjj > 1.4 PeV, CiT <0.8 and rfgtr < 60 m gives 
the best correlation with the Galactic Plane. Because of the 
small size of the burn sample, this correlation is most likely 
accidental. 

The optimized cuts are applied to the complete IC40 data 
set minus the burn sample. There are 268 candidate events 
of which 28 are located in the source region. Fig. 2 is 
a map of the sky showing all 268 events. The colors in 
the background are the integrated neutral atomic hydro- 
gen (HI) column densities based on data from the Lei- 
den/Argentine/Bonn survey [6]. These values have no im- 
pact on the analysis and are only plotted to indicate the 
positions of the Galactic Plane and the Magellanic Clouds. 
The significance of the correlation with the Galactic Plane 
is again tested by producing data sets with scrambled RA. 
An equal or higher number of source region events is found 
in 2 1 % of the scrambled data sets . The source region has a 
non-significant upward fluctuation of +0.9a. 

We follow the procedure of Feldman and Cousins [8] to 
construct an upper limit on the ratio of gamma rays to CRs. 
The background is determined by selecting a range of RA 
that avoids the source region. Within this range the data 
points are scrambled multiple times and for each scrambled 
set the number of events in a pre-defined region of the same 
shape and size as the source region is counted. This yields 
a mean background of 24.13 events for the source region. 
To derive at an upper limit on the ratio between the flux of 
gamma rays to CR several quantities have to be included 
which have been found with MC studies. The energy re- 
construction of IceTop is cahbrated on proton showers and 
overestimates the gamma-ray shower energy by 14%. The 
lower energy limit of 1 .4 PeV on CR showers corresponds 
to 1 .2 1 PeV gamma-ray showers . Other quantities that have 
been found with simulation are the fraction of gamma-ray 
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Figure 3: Limits on the gamma-ray to CR flux ratio from 
wifliin 10 degrees of the Galactic Plane by CASA-MIA [3] 
and 1C40 analysis. The sensitivity of 5 years fuU IceCube 
data is indicated by the dashed line. 

showers producing a muon signal in IceCube, and the dif- 
ference in the effective area for CR and gamma ray show- 
ers. Finally we arrive at a 90% C.L. upper hmit on the ra- 
tio of gamma-to-CR shower above 1.21 PeV of 8.1 • 10~^, 
plotted in Fig. 3. 

32 Unbiimed point source search 

In an additional search for point-like sources we account 
for the possibility that a single source dominates the PeV 
gamma ray sky and that this source does not necessarily 
lie close to the Galactic Plane. An unbinned point source 
search is performed on the sky within the declination range 
of -80° to -60° in a method that foUows [9]. 

The data is described by an unknown amount of signal 
events on top of a flat distribution of background events. 
In an unbinned search, a dense grid of points in the sky is 
scanned. For each point a maximum UkeUhood fit is done 
for the relative contribution of source events over back- 
ground events. 
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Figure 4: Lambda contours of an unbinned search for gamma-ray point sources. 



For a particular event i in a sample of N events, the proba- 
bihty density function (PDF) is given by 

Pdns) = ^S,+ [l-^)B, (1) 

where ns is the number of events that is associated to the 
source, Bi is a flat background PDF, and 

is the two-dimensional Gaussian source PDF, in which A^" 
is the space angle between the event and the source test 
location, and cr = 1.5° is the angular resolution of IceTop. 
For each point in the sky there is a Ukehhood function 

L{ns)=UiPi{ns), (3) 

and associated test statistic 

A = -2 (log(L(0) - L{ns)) , (4) 

which is optimized for ng. 

Fig. 4 is a map of the sky within range of -80 ° to -60° show- 
ing the events in this region and contours of the test statis- 
tic A. The maximum value is A = 5.5 at 6 = —1.34 and 
RA= 3.99, corresponding to 7.25 signal events. The sig- 
nificance of this value for A is found by producing 10, 000 
scrambled data sets by randomizing the RA of each data 
point. The mean test statistic value for the hottest spot is 
A = 5.74, so the actual data set has a shght underfluctua- 
tion of -O.lcr and is fully consistent with a flat background. 
A 3a upward fluctuation would have corresponded with 
A = 12.5. 

An upper limit on the gamma-ray flux from point sources 
can be set by assuming all 7.25 events associated to the 
hottest spot are actual gamma-rays from a point source. 
The effective area of IC40 for a source at a zenith angle 
of 13.2° is 5.7 • 10^ m^ and the total observation time is 
3.2 • 10*^ s. Assuming an spectrum with index 7 = —2 an 
upper limit of 7.0 • IQ-^^ cm-^s-^TeV-^ at £ = 1.2 PeV 
can be placed with a C .L. of 99% . 

4 Full IceCube sensitivity 

The construction of IceCube was finished in the winter of 
20 10/11. With a larger array the acceptance increases con- 



siderably. Because of the condition that the shower axis 
has to be inside the instrumented area of both IceCube 
and IceTop, the increase is especially dramatic at higher 
zenith angles. The maximum zenith angle is also larger, 
which extends the FOV to cover a larger part of the Galac- 
tic Plane and probe an area even closer to the Galactic 
Center. Additional efficiency is gained by using the in- 
formation of isolated hits. In the IC40 configuration, hit 
information was only available when neighboring DOMs 
both detected a signal, while isolated hits were discarded. 
The isolated hit information will increase the discriminat- 
ing power between CRs and gamma-rays. Furthermore, 
progress in shower core reconstruction algorithms have re- 
duced the number of mis-reconstructed showers. The ex- 
pected sensitivity of 5 years of IceCube data is plotted in 
Fig. 3. 

5 Conclusion 

We have investigated the sensitivity of IceCube for the de- 
tection of PeV gamma-rays from the Galactic plane. With 
one year of data from the unfinished detector we have de- 
rived a limit on the flux of gamma rays to CR which is not 
competitive with existing limits. The full IceCube detector 
wiU be able to reach a sensitivity to the diffuse Galactic flux 
similar to that of C ASA-MI A, for slightly higher energies. 
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Abstract: We report on a study of the energy dependence in the arrival direction distribution of cosmic rays at median 
energies of 20 and 400 TeV. The data used in this analysis contain 33 x 10^ downward going muon events collected 
by the IceCube neutrino observatory between May 2009 and May 2010 when it comprised of 59 strings. The high 
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1 Introduction 

During the last decades. Galactic cosmic rays have been 
found to have a small but measurable energy dependent 
sidereal anisotropy in their arrival direction distribution 
with a relative amplitude of order of 10"'' to 10"*^. The 
first comprehensive observation of the cosmic ray side- 
real anisotropy was provided by a network of muon detec- 
tors sensitive to cosmic rays between 10 and several hun- 
dred GeV [1]. More recent underground and surface ar- 
ray experiments in the Northern hemisphere have shown 
that a sidereal anisotropy is present in the TeV energy 
range (i.e. Tibet Air Shower Arrays (Tibet AS7) [2], Mi- 
lagro [3]). The IceCube neutrino observatory reported for 
the first time the observation of the cosmic rays sidereal 
anisotropy in the Southern sky at energies in excess of 
about 10 TeV [4]. The cosmic ray anisotropy reported by 
IceCube showed that the large scale features were a contin- 
uation of those observed in the Northern hemisphere in the 
same energy range. 

At high energies, the Tibet AS7 collaboration reported 
the non-observation of a sidereal anisotropy in the cosmic 
ray arrival direction distribution at ^ 300 TeV [2]. An- 
other measurement at ~370 TeV was recently published 
by the EAS-TOP collaboration reporting a larger sidereal 
anisotropy in amplitude with a shift in phase from that ob- 
served at lower energies [5]. At this point the two obser- 



vations do not provide a coherent picture of the sidereal 
anisotropy at high energy in the Northern hemisphere. 

In this contribution we present the analysis and the results 
of the study of large scale cosmic ray anisotropy at median 
energies of 20 and 400 TeV by IceCube for the Southern 
celestial sky. Observations of the energy dependence of the 
anisotropy could provide us with further information for 
the development of theoretical models to better understand 
the origin and propagation of cosmic rays. 

2 Data Analysis 

The main goal of the IceCube neutrino observatory is to 
detect astrophysical neutrinos. Neutrinos passing through 
the Earth and interacting in the vicinity of IceCube produce 
muons or other secondaries that emit Cherenkov radiation 
in the clear ice surrounding the detector. It is these light 
signals that the IceCube optical modules record. On the 
other hand, muons produced by cosmic ray air-showers in 
the atmosphere above IceCube are also detected. These are 
observed as down-going tracks. 

The data used in this analysis are the downward going 
muons collected by the IceCube detector comprising 59 
strings. The data were collected from May 2009 to May 
2010 [6]. The events used in this analysis are those re- 
constructed by an online likelihood based reconstruction 
algorithm at an average rate of -^1250 Hz with a median 
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angular resolution of ^ 3 ° . A range of selection criteria is 
applied to that data to ensure good quality and stable runs. 
The final data set consists of 33 x 10^ events corresponding 
to a detector livetime of 324.8 days. 

Due to the detector's unique location, the IceCube Obser- 
vatory has full coverage of the Southern sky up to declina- 
tion angles of -25° degrees at any time of the year. The sea- 
sonal and atmospheric variation of the muon intensity were 
found to occur uniformly across the entire field of view, 
and therefore it did not affect the arrival distribution of the 
reconstructed events [4]. The crucial effect that needed to 
be accounted for was due to the geometrical shape of Ice- 
Cube: events that were parallel to a larger number of strings 
were more efficiently reconstructed. The non-uniform time 
coverage due to detector downtime and run selection pre- 
cluded the complete averaging of the detector geometrical 
asymmetry and generated artificial variations of the arrival 
direction of cosmic rays in equatorial coordinates. The lo- 
cal azimuthal asymmetry was corrected by reweighting the 
arrival directions of the data [4]. 

3 Energy Estimation 

The energy dependence of the large scale anisotropy may 
hint at the nature of the source (or sources) of the cosmic 
rays, as well as their propagation through galactic magnetic 
fields. Similar measurements have been carried out over a 
wide range of energies by underground muon experiments 
and air shower arrays [1,7,8,9, 10] but there are relatively 
few observations in the multi-TeV region. 

Since IceCube detects cosmic ray properties indirectly 
through the observation of muons produced in the exten- 
sive air showers, the cosmic ray particles energy is inferred 
from the estimation of the muon energy. In this analysis, 
we use the number of optical modules Nch participating 
together with the zenith angle 9 of the event to estimate the 
energy of the events. Using Monte Carlo simulations, we 
have identified cuts of constant energy in (Nch, 6) for the 
two event samples. The primary cosmic ray energy estima- 
tion had a resolution of about 0.5 in log scale and, this is 
due to the fact that it is dominated by the large fluctuations 
of the number and energy of muons produced in the exten- 
sive air showers. The low energy sample contained events 
with a median energy of 20 TeV, where 68% of the events 
were between 4 — 63 TeV; and the high energy sample con- 
tained events with a median energy of 400 TeV, were 68% 
of the events were between 100 — 1258 TeV. 

4 Results 

To investigate the arrival direction distribution of the cos- 
mic rays, we study the relative intensity of the cosmic ray 
flux. The arrival direction distribution is dominated by the 
zenith angle dependence of the flux. Therefore, we normal- 
ize the flux within declination belts of width ^ 3 ° , which 
corresponds to the angular resolution of the data. 



To quantify the scale of the anisotropy, we fitted the right 
ascension projection distribution for declinations angles 
-25°to -75° degrees of the data to a first and second-order 
harmonic function of the form 

2 

Y,AjCos[j{a-<i>j)]+B (1) 
i=i 

where [Aj , 0i ) are the amplitude and phase of the 
anisotropy, a is the right ascension, and _B is a constant. 
Figures 1 and 2 show the relative intensity cosmic ray 
maps together with the profile of the data in right ascen- 
sion for events in both the 20 TeV and the 400 TeV en- 
ergy sample. The error bars in the right ascension projec- 
tion are derived by propagating the statistical errors from 
each declination belt, and the gray band indicates the esti- 
mated maximal systematic uncertainties (described in sec- 
tion 4.1.1). The solid line indicates the fit of eq. (1) to the 
data. The first and second harmonic amplitude and phase of 
the sidereal anisotropy for the low energy sample together 
with their statistical and systematic uncertainties are Ai = 
(7.9±O.l,tot.±O.4,j,,t.)xlO-4and0i = 50.5°±1.0°,„,± 
l.l°y,t,,A2 = (2.9 ± OAstat. ± OA,yst.) X 10-^ and 02 = 
299.5° ± 1.3°tat ± ^-^syst.- While those for the high en- 
ergy sample are A i = (3.7±0.7stat ±0.7 syst )xlO~**and 
= 239.2° ±10.6°,„,± 10.8°^,,, A2 = {2.7 ±0.7, tat ± 
0.6,,,O X 10-4 and <^2 = 152.7° ± 7.0°,„, ± 4.2°^,,,. 

4.1 Reliability Checks 
4.1.1 Data Stability 

In order to assess and quantify the systematic uncertainties 
in the sidereal anisotropy of cosmic ray arrival direction 
distribution a number of checks were applied by dividing 
the low and high energy data samples in exclusive halves 
based on different criteria. A full analysis was then ap- 
plied for each dataset and the relative intensity distribution 
in right ascension was determined for each of them. The 
stability checks applied are: 

• Seasonal variations dependence: where the data was 
divided in two seasons (winter and summer). 

• Rate variations dependence: where the data was di- 
vided for each day by the sub-run (a sub-run corre- 
sponded to approximately 2 minutes of observations) 
rate fluctuating being greater than or less than the 
median rate of the day. 

• Choices of events sample dependence: where the 
data was divided by the sub-run number using multi- 
ple categories. 

• Non-uniform time coverage dependence due to de- 
tector down time and quality run selection: where the 
analysis was performed on the sub-sample of days 
with maximal data collection time. 
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The sidereal distribution of relative intensity in the arrival 
direction of the cosmic rays for the low and high energy 
samples are used to evaluate the spread in the experimen- 
tal observation from the full-year event samples. The gray 
bands in the right ascension projection in Figures 1 and 2 
describe the maximal spread obtained from the result of all 
the stability checks described in this section. 

4.1.2 Solar Diurnal Anisotropy 

To test for the stability of the observatory and its time cov- 
erage, an effective way to have an absolute calibration of 
the experimental sensitivity for the detection of the sidereal 
directional asymmetries is to measure the solar anisotropy 
from the Earth's revolution around the Sun. This observa- 
tion has solid theoretical grounds and it was first reported 
in 1986 [11] and then later observed by multiple experi- 
ments in the multi-TeV energy range (i.e. [2], [3]). The 
observed solar anisotropy is consistent with a dipole that 
describes an apparent excess of cosmic rays in the direc- 
tion of Earth's motion around the Sun and a deficit in the 
opposite direction. 

Figure 3 shows the projections of the cosmic ray arrival 
direction in solar reference frame, for both energy sam- 
ples (20 and 400 TeV). The error bars are the statistical 
errors, and the shaded bands indicate the expectation of the 
dipole. A fit to the projection of relative intensity distribu- 
tion is shown by the black line and was done using the first 
harmonic term of eq. (1). The figures show that the experi- 
mental observation of the solar dipole is consistent with the 
expectations in both amplitude and phase. This observation 
of the solar diurnal anisotropy supports the reliability of the 
sidereal anisotropy determination. 

4.1.3 Anti-sidereal Anisotropy 

An annual modulation in the amplitude of the solar 
anisotropy is expected to result in a spurious effect in the 
sidereal anisotropy. This would produce a bias in the ob- 
served sidereal anisotropy. To estimate this bias the so- 
called anti-sidereal time, i.e. a non-physical time frame ob- 
tained by reversing the sign of the transformation from so- 
lar time to sidereal time [12] . 

Figure 4 shows the projections of the cosmic ray arrival 
direction in anti-sidereal reference frame, for both energy 
samples (20 and 400 TeV). The error bars are the statisti- 
cal errors. The distributions were then fitted to the dipole 
term of eq. (1). The uncertainty in the first harmonic am- 
plitude and phase implied by the study in the anti-sidereal 
time frame is within the statistical and systematic errors de- 
termined from the data stability tests. The absence of the 
signal in the anti-sidereal time insures the reliability of the 
anisotropy observed in sidereal time. 
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Figure 1: The top figure is the IceCube relative intensity 
cosmic ray map for the low energy sample (Median energy 
of the primary cosmic ray particle of 20 TeV). The bottom 
figure is the one dimensional projection in right ascension 
a of the two-dimensional cosmic ray map. The black line 
corresponds to the first and second harmonic fit to the data. 
The gray band indicates the estimated maximal systematic 
uncertainties. 
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Figure 2: The top figure is the IceCube relative intensity 
cosmic ray map for the high energy sample (Median energy 
of the primary cosmic ray particle of 400 TeV) . The bottom 
figures is the one dimensional projection in right ascension 
a of the two-dimensional cosmic ray map. The black line 
corresponds to the first and second harmonic fit to the data. 
The gray band indicates the estimated maximal systematic 
uncertainties. 

5 Conclusion 



In this contribution we presented the results on the large 
scale cosmic ray sidereal anisotropy at cosmic ray median 
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Figure 3: The solar projections for cosmic rays with me- 
dian energy of 20 TeV (The top figure) and 400 TeV (The 
bottom figure). The error bars are the statistical errors, and 
the shaded band indicated the solar diurnal dipole expected 
from the motion of the Earth around the Sun. 



Figure 4: The anti-sidereal projections for cosmic rays with 
median energy of 20 TeV (The top figure) and 400 TeV 
(The bottom figure). The error bars are the statistical errors. 
The black line corresponds to the first harmonic fit to the 
data. 



energies of 20 TeV and 400 TeV. The source of the sidereal 
anisotropy is still unknown. It is believed that a possible 
contribution to this observed anisotropy might be from the 
Compton-Getting effect, due to the orbital motion of the so- 
lar system around the galactic disk. However, the sidereal 
anisotropy from both energy samples do not appear to be 
consistent with that expected from the suggested Compton- 
Getting model neither in amplitude nor in phase. 

The sidereal anisotropy observed at 20 TeV with IceCube- 
59 is consistent with the previously reported observation 
with IceCube [4] , thus providing a confirmation of an ap- 
parent continuation of the arrival distribution pattern ob- 
served in the Northern hemisphere. On the other hand the 
sidereal anisotropy observed at 400 TeV shows substan- 
tial differences with respect to that at lower energy. The 
anisotropy at high energy shows a relative deficit region 
in right ascension where the broad excess dominated at 
primary median energy of 20 TeV. Also the wide relative 
deficit region at low energy seemed to have disappeared 
at primary median energy of 400 TeV. Whatever generated 
the sidereal anisotropy at 20 TeV seems to have no effect at 
400 TeV. 

This is the first observation of the sidereal anisotropy at 400 
TeV in the Southern hemisphere. We are continuously an- 
alyzing events from IceCube with updated configurations. 
IceCube construction is now completed with 86 strings de- 
ployed with a volume of km'^ in January of 201 1 . With the 



higher statistical power expected from the observed cosmic 
ray muons we will be able to improve our understanding of 
the energy dependence of the anisotropy with more signifi- 
cance and close to the knee region. The energy dependence 
of the cosmic ray anisotropy is vital to our understanding 
of the source and propagation of cosmic rays. 
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1 Introduction 

The IceCube detector, deployed between 1450 m to 2450 m 
below the surface of the South Polar ice sheet, is de- 
signed to detect upward-going neutrinos from astrophysical 
sources. However, it is also sensitive to downward-going 
muons from cosmic-ray air showers. To penetrate the ice 
and trigger the detector, the muons must be produced by 
cosmic rays with energies of at least several TeV. The ar- 
rival direction of a cosmic muon is typically within 0.2° of 
the direction of the primary particle, so the arrival direction 
distribution of muons is also a map of cosmic ray arrival 
directions between about 1 and several hundred TeV. 

At energies of a few TeV, it is believed that galactic 
magnetic fields should randomize the arrival directions 
of charged cosmic rays. However, in recent years an 
anisotropy in the arrival direction distribution has been 
reported on both large and small angular scales. The 
anisotropy is observed in the northern sky between sev- 
eral TeV and several hundred TeV by the Tibet AS 7 array 
[1], Super- Kamiokande [2], Milagro [3, 4], ARGO-YBJ 
[5], and EAS-TOP [6]. In 2010, an analysis of muons 
recorded by the IceCube detector also revealed a large- 
scale anisotropy in the southern sky [7]. In these proceed- 
ings, we present the results of a search of the southem sky 
for anisotropy on all angular scales using data recorded be- 
tween May 2009 and May 2010. 



2 The IceCube Detector and Data 

IceCube is a km'^-size neutrino detector frozen into the 
glacial ice sheet at the geographic South Pole. High- 
energy charged particles passing through the detector emit 
Cherenkov radiation, and their tracks are recorded by an 
array of Digital Optical Modules (DOMs) embedded in the 
ice. The DOMs are attached to 86 vertical cables, or strings 
(with 60 DOMs per string), which are used to transmit data 
to the surface. Construction of IceCube was completed in 
December 2010. The detector has been operating in vari- 
ous configurations since construction began. Between 2007 
and 2008, it was operated with 22 strings deployed (IC22), 
between 2008 and 2009 with 40 strings (IC40), and be- 
tween 2009 and 2010 with 59 strings. 

Muons are identified using a simple majority trigger, which 
requires 8 or more DOMs in local coincidence within a 5 /iS 
window. The trigger rate of downgoing muons is 0.5 kHz 
in IC22, 1.2 kHz in IC40, and 1.4 kHz in IC59, about a fac- 
tor of 10^ larger than the neutrino rate. The rate is too large 
to transfer complete waveforms for all events via satellite, 
so the muon data are reconstructed on-line and compressed 
using a Data Storage and Transfer (DST) format. Muon 
tracks are identified using a maximum-likelihood recon- 
struction, and the event arrival direction and energy estima- 
tors are stored in DST files and sent north by satellite. Dig- 
itized waveforms are discarded due to limited bandwidth. 
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The analysis presented in this paper uses the DST data col- 
lected during IC59 operations between May 20, 2009 and 
May 30, 2010. The data set contains nearly 3.4 x 10^° 
muon events recorded during a live time of 335.5 days. A 
cut in zenith angle at 65° removes less well-reconstructed 
tracks near the horizon, reducing the final analysis sample 
to 3.2 X 10^° events. The median angular resolution of 
events in the sample is 3°; unlike the neutrino analysis, the 
cuts have not been optimized for point source identifica- 
tion. Using simulated data, we estimate the median energy 
of the events to be 20 TeV. The energy resolution of the 
cosmic ray primaries is about A log (E/eW) — 0.5. 

3 Analysis 

3.1 Calculation of the Reference Level 

The arrival direction distribution of cosmic rays in IceCube 
exhibits anisotropy caused by non-physical effects such as 
gaps in the detector uptime. These effects must be re- 
moved before any physical anisotropy can be identified. 
Therefore, the first step in the analysis is the creation of 
a "reference map," a sky map that describes what the ar- 
rival direction distribution would be if the cosmic ray flux 
were isotropic. The reference map must be subtracted from 
the arrival direction distribution in the data to find regions 
where the cosmic ray flux deviates from the isotropic ex- 
pectation. 

We estimate a reference map for IC59 data using the time- 
scrambling method of Alexandreas et al. [8]. The sky 
is binned into an equal-area grid in equatorial coordinates 
with 0.9° resolution using the publicly-available HEALPix 
library [9]. From this sample two sky maps are produced. 
The data map N{a, S) stores the arrival directions of all 
events. For each detected event in the data map, "fake" 
events are generated by keeping the local zenith and az- 
imuth angles {6, (f)) fixed but reassigning the time to that of 
another event recorded within a pre-defined time window 
At. The event times within the window are resampled 20 
times and used to generate 20 fake celestial coordinate pairs 
{a, S) which are used to fill a reference sky map {N{a, S)) 
with a weight of 1/20. 

A known disadvantage of the method is that it can be bi- 
ased by a strong anisotropy, leading to artificial deficits or 
excesses next to regions of true excess or deficit [3]. How- 
ever, the procedure does preserve the local arrival direction 
distribution of the data and naturally compensates for vari- 
ations in the event rate which are difficult to model, such as 
changes in atmospheric conditions which affect the produc- 
tion of muons in air showers. The only critical requirement 
for time scrambling is to choose the buffer length At such 
that detector conditions remain stable during this period. 
Using tests to compare the distribution of arrival direc- 
tions in local coordinates across measurement periods, we 
have verified that the IceCube detector is stable over pe- 
riods of at least one full day. In this analysis we choose 
At = 24 hr. 
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Figure 1 : Top: Relative intensity of the IC59 data produced 
with At = 24 hr. Middle: Fit of dipole and quadrupole 
moments to the relative intensity. Bottom: Residual map 
after removal of the dipole and quadrupole. 

3.2 Relative Intensity 

Once the data and reference maps are calculated, deviations 
from isotropy can be analyzed by computing the relative 
intensity 

AN, ^ N,ia,d)-{N,ia,S)) 

(TV),; {N,{a,S)) • ^ ^ 

which gives the amplitude of deviations from the isotropic 
expectation in each angular bin i. The significance of the 
deviation in bin i can be calculated using the method of Li 
and Ma [10]. 

A map of relative intensity of the IC59 data, binned with 
0.9° resolution, is shown at the top of Fig. 1. The map ex- 
hibits obvious correlations between bins, such as a broad 
excess in the relative counts near a = 105° and a broad 
deficit near a = 225°. The relative intensity in these re- 
gions is of order 10~^. This structure is the large-scale 
anisotropy reported in the IC22 data by Abbasi et al. [7]. 
Since the IC59 data set is larger than that of IC22 by an or- 
der of magnitude, it is possible to see the large-scale struc- 
ture in the data without further rebinning. 

3.3 Removal of Large-Scale Structure 

The relative intensity shown at the top of Fig. 1 is dom- 
inated by large-scale structures, but there are additional 
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Figure 2: IC59 residual maps after subtraction of dipole and quadrupole moments, showing relative intensity (left) and 
significance before trial factors are applied (right). The maps have been smoothed with a circular window of 20°. 
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Figure 3: Residual maps from IC22 (left) and IC40 (right) showing the significance of deficit and excess regions with 20° 
smoothing. 



small-scale features in the map. This can be demonstrated 
by calculating the angular power spectrum of the data, 
which exhibits significant power not only at the largest an- 
gular scales, but also down to scales of about 10° (see de- 
tailed discussion in Abbasi et al. [11]). 

We can also demonstrate the presence of small-scale struc- 
ture by explicitly removing the largest angular correlations 
in the data. This is accomplished by fitting dipole and 
quadrupole terms to the relative intensity map and then 
subtracting the fit to obtain a map of residual counts. In 
the middle panel of Fig. 1, we show the result of the fit 
of dipole and quadrupole moments to the intensity map. 
By themselves, these two terms account for much of the 
amplitude of the per-mille anisotropy observed in the data. 
However, the fit x^/ndf = 14743/14187 corresponds to a 
X^-probability of approximately 0.05%, suggesting that the 
dipole and quadrupole are not sufficient to explain all of the 
structures observed in the angular distribution of AN/ (N) . 

Subtraction of the dipole and quadrupole fit from the data 
gives the bottom panel in Fig. 1. The fit residuals are rel- 
atively featureless at first glance, but the bin size is not 
optimized for a study of significant anisotropy at angular 
scales larger than the angular resolution of the detector To 
increase the sensitivity to the small-scale structure in the 
data, we apply a smoothing procedure which takes the ref- 
erence level and residual data counts in each bin and adds 
the counts from pixels within some angular radius of the 
bin. This procedure results in a map with Poisson uncer- 
tainties, though the bins are no longer statistically indepen- 
dent. 

Smoothed residual maps of relative intensity and signifi- 
cance are shown in Fig. 2. To make this figure, a smooth- 



ing radius of 20° was chosen. Strong regions of excess 
and deficit are visible in the data. Note that the actual size 
of any given excess or deficit region (and thus the optimal 
smoothing scale) is not known a priori. Therefore, we study 
the sky map on all smoothing scales from 3° (the angular 
resolution of the data) to 45° in steps of 1° and search for 
regions of high significance at any location. 

Applying this procedure, the two most significant features 
on the sky are a region with a peak significance of 5.3cr at 
a smoothing radius of 22° {a = 122.4°, (5 = -47.4°) and 
a region of peak significance 4.9cr at a smoothing radius 
of 13° (a = 263.0°, -44.1°). The significance values ac- 
count for statistical trials due to the scan over smoothing 
radii and the location of the most significant pixels. The 
trial factors were estimated using a Monte Carlo simula- 
tion of an isotropic flux recorded by the IceCube detector. 

4 Systematic Checks 

Several tests have been performed on the data to ensure the 
stability of the observed anisotropy and to rule out possi- 
ble sources of systematic bias. Among the influences that 
might cause spurious non-physical anisotropies are the de- 
tector geometry, the detector livetime, non-uniform expo- 
sure of the detector to different regions of the sky, and diur- 
nal and seasonal variations in atmospheric conditions. Due 
to its unique location at the South Pole, many of these ef- 
fects play a lesser role for IceCube than for detectors at the 
middle latitudes. We check the validity of this assumption 
by searching the data for a signal expected in solar coordi- 
nates, and by examining the data from previous configura- 
tions of the detector. 
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4.1 Solar Dipole Analysis 

The cosmic rays at TeV energies do not co-rotate with the 
Earth about the Sun, and so it is expected that the flux of 
cosmic rays should exhibit a dipole modulation in solar 
coordinates. The expected change in relative intensity is 
given by 

^ = (7 + 2)-cos^^, (2) 
I c 

where / is the cosmic ray intensity, 7 = 2.7 is the power 
law index of the energy spectrum at several TeV, v/cis the 
ratio of the orbital velocity of the Earth with respect to the 
speed of light, and is the angle between the cosmic ray 
arrival direction and the direction of motion of the Earth. 
Given v = 30 km s~^, we expect an amplitude of 4.7 x 
10-4. 

The solar dipole provides a convenient measurement with 
which to check the analysis technique. The reference 
and data maps were computed in solar coordinates, and 
the relative intensity map has been fit with a dipole func- 
tion. A single dipole describes the data quite well: the fit 
xVndf = 14206.8/14192 corresponds to a x^-probabiUty 
of 41.6%. The dipole is aligned at a longitude of 270° in 
solar coordinates, following the expectation for a relative- 
motion anisotropy in this coordinate frame. Its amphtude 
is (3.66 ± 0.14 ± 0.99) x 10"^, in agreement with flie ex- 
pectation within the large systematic uncertainty. A more 
detailed study of the solar dipole anisotropy is presented in 
Abbasi and Desiati [12]. 

4.2 Anti-Sidereal Time Analysis 

Seasonal variations in the solar dipole anisotropy can create 
spurious signals in equatorial coordinates, and vice-versa. 
We study this effect using an artificial time scale called 
anti-sidereal time, which is calculated by reversing the sign 
of the transformation between universal time and sidereal 
time. No physical anisotropy is expected in this reference 
frame, but it can be used to identify systematic distortions 
due to seasonal effects. 

We have produced sky maps using anti-sidereal coordinates 
and performed the same analysis presented in Section 3. 
No regions of significant excess or deficit are observed in 
the anti-sidereal maps for any smoothing scale. A detailed 
discussion is provided in [1 1]. 

4.3 Comparison to IC22 and IC40 Data 

An important cross-check of the 1C59 analysis can be made 
by applying the same method to data recorded in the two 
data periods prior to IC59. The IC22 data set contains 
4 X 10^ events recorded between July 2007 and April 2008, 
and the IC40 data set contains 1.9 x 10^° events recorded 
between April 2008 and April 2009. While the samples are 
smaller than IC59, we expect to observe the most promi- 
nent structures in these data, albeit with reduced signifi- 
cance. The IC22 and IC40 data can be used to verify that 



the structures observed in the arrival direction distribution 
do not depend on the geometry of the detector or the data- 
taking period. For example, the shape of 1C22 was highly 
asymmetrical, with a long and a short axis. 

In Fig. 3, we show the residual maps from 1C22 and 1C40 
after the subtraction of dipole and quadrupole terms and 
smoothing by 20°. While none of the features in IC22 and 
IC40 have a pre-trial significance above 5a, they align with 
the regions of deficit and excess observed with IC59 data 
(cf. Fig 2). The main features on both large and smaU 
scales appear to be persistent in all data sets. 

5 Conclusions 

Using 3.2 X 10^° events recorded with the partially- 
deployed IceCube detector between May 2009 and May 
2010, we have found that the arrival direction distribu- 
tion of cosmic rays at several TeV exhibits significant 
anisotropy on several angular scales. The data are domi- 
nated by dipole and quadrupole moments, but there is also 
significant structure on angular scales down to about 10°. 
These structures become visible in the sky map when the 
dominant terms are subtracted. 

There is currently no explanation for these local enhance- 
ments in the cosmic ray flux, and so the study of cosmic 
ray arrival directions in the TeV region will continue to be a 
major ongoing effort in IceCube. We are currently studying 
the anisotropy in the 79-string configuration of the detector 
(IC79). During the next several years, with the IceCube de- 
tector operating in its complete 86-string configuration, the 
data will increase by 4.5 x 10^° events per year. With these 
high statistics we will be able to study possible time depen- 
dencies in the anisotropy and compare to similar studies in 
the northern hemisphere [4, 13] 
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Abstract: 

IceCube is a kilometer scale neutrino observatory that collects a large number of cosmic ray induced muon events. These 
events, which are background for neutrino searches, are observed at a rate that is suitable for high-statistics studies of 
cosmic rays in the Southern hemisphere. The apparent anisotropy created by the motion of the Earth around the Sun, the 
solar dipole, is systematically analyzed. The solar dipole is simulated, and the predictions for the integrated effect over an 
entire year and over shorter periods of a quarter-year are compared to data. The experimental observation is found to be 
in good agreement with the expectation. Finally, we show that the interference between the solar dipole and the sidereal 
anisotropy is well understood within the statistical uncertainties. 
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1 Introduction 

The motion of an observer in the cosmic ray plasma is pre- 
dicted to cause an apparent dipole effect in the cosmic ray 
arrival direction; this effect was first pointed out in 1935 by 
Arthur Compton and Ivan Getting [1]. Using this principle, 
the motion of the Earth around the Sun would produce an 
excess in the direction of motion and a deficit in the oppo- 
site direction. This observation was first reported by Cutler 
in 1986 [2], and then later observed by multiple experi- 
ments (i.e. Tibet Air Shower Arrays (Tibet AS7) [3], Mila- 
gro [4], and EAS-TOP [5]). This measurement is referred 
to in this contribution as the solar dipole effect anisotropy 
and is expressed as 

AJ V 

— = (7 + 2)-cos(0,) (1) 
(7) c 

where / is the intensity, 7 is the differential cosmic ray 
spectral index, v is the Earth's velocity, c is the speed of 
Ught, and 6^ is the angle between the reconstructed arrival 
direction of the cosmic rays and the direction of the motion 
of the Earth [6]. 

In addition to the solar dipole effect, galactic cosmic ray 
anisotropies have been observed by underground and sur- 
face array experiments in the Northern hemisphere such 
as Tibet AS7 [3]. Furthermore, IceCube recently reported 
the first observation of the sidereal anisotropy in cosmic 
ray arrival direction in the Southern sky [7] with a relative 



ampUtude of the order of 10~^ to 10~^. The cosmic ray 
anisotropy showed that the sidereal anisotropy in the South- 
ern sky is a continuation of that observed in the Northern 
hemisphere. 

In this proceeding, we report on the study of the observa- 
tion and the expectation of the solar dipole effect in one 
full year. We also examine the effect of the mutual interfer- 
ence between the sidereal anisotropy and the solar dipole 
effect by examining quarter-year periods. The solar dipole 
effect and the interference between the sidereal and the so- 
lar effect over quarter-year periods is found to be well un- 
derstood. This understanding supports the reliability of the 
sidereal anisotropy observations. 

2 Data Analysis 

IceCube is optimized for the detection of up-going, high- 
energy neutrinos, however its trigger rate is dominated by 
down-going atmospheric muons created in cosmic ray air- 
showers above the detector. In this work, we use the down- 
going atmospheric muon flux as our signal to study cosmic 
rays in the Southern sky. The data used in this analysis are 
the down-going muons collected by the IceCube detector 
in its 59-string configuration from May 2009 to May 2010. 
The data are reconstructed by an online maximum likeli- 
hood based reconstruction algorithm at an average rate of 
~1.2 kHz. A range of selection criteria is appUed to the 
data to ensure detector stabihty. The final data set consists 
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Figure 2: The plot shows the one-dimensional projection in 
Q! — ctsun of the relative intensity of cosmic rays. The data 
are shown with statistical uncertainties in black crosses, 
and the black line corresponds to the expectation from 
eq. 1. The gray shaded band indicates the 68% spread of 
the uncertainty in the simulated solar dipole. 



of 33 billion events and corresponds to a detector livetime 
of 324.8 days. The events are initiated by muons from cos- 
mic ray air-showers with a median primary energy of 20 
TeV. The angular resolution between the muon to the pri- 
mary cosmic ray particle has a median value of ~3°. 

Located at the South Pole, IceCube observes the Southern 
sky year-round. However, an artificial azimuthal asymme- 
try occurs as an artifact of the detector's geometrical con- 
figuration (Figure 1) combined with non-uniform time cov- 
erage. This azimuthal asymmetry was eliminated by re- 
weighting the azimuthal arrival directions of the incoming 
cosmic rays [7]. 

3 Solar Dipole Anisotropy Measurement and 
Expectation 

The solar dipole is a well-known signal, both from the- 
ory and from observations by other experiments. A mea- 
surement of this well-known signal with IceCube strongly 
supports the reliability of the observation of the sidereal 
anisotropy with the same detector [8]. 

The amplitude of the solar dipole depends on the geograph- 
ical latitude of the observer and on the angular distribution 
of the detected cosmic ray events at the observatory. Due to 
the location of the IceCube Observatory at the South Pole, 
the sky is fully visible at any given time. Therefore, the 
solar dipole is observed in a reference system where the lo- 
cation of the Sun is fixed, where the latitude coordinate is 
the declination, and the longitude coordinate is defined to 
be the difference between the right ascension of the cos- 
mic ray arrival direction and the right ascension of the Sun 
(a — asun)- In this reference frame, the excess due to the 
solar dipole is expected to be at 270° and a deficit at 90°. 

To verify that the experimental observation of the solar 
dipole is consistent with the expectation, the predicted pro- 
jection of the solar anisotropy was calculated for the Ice- 



Cube location. Since the observed amplitude of the solar 
dipole depends on the actual angular distribution of the de- 
tected muon events, the calculation of the expectation was 
performed using the azimuth-corrected experimental data. 
The expectation of the solar dipole was calculated by com- 
puting the relative intensity of the solar dipole through the 
cosmic ray plasma (through eq. 1). Instead of counting the 
number of events within a given bin in right ascension from 
the Sun, we calculated for each event, a mean weight cor- 
responding to the expected relative intensity of the solar 
dipole. 

The uncertainties in the cosmic ray spectral index, in the 
Earth's velocity, and in the reconstructed arrival direction 
of the events were included into the calculation of the 
uncertainty of the expectation. The mean spectral index 
was evaluated using the all-particle cosmic ray spectrum 
from [9] and was found to be 2.67 ± 0.19. The value used 
for Earth's velocity was v = 29.8 ± 0.5 km/s, where the 
error takes into account the spread between the maximum 
and minimum along the elliptical orbit. The angle 6y be- 
tween the reconstructed direction of the muon events and 
the Earth's velocity vector at the time the event was de- 
tected was evaluated accounting for the experimental point 
spread function. The expected solar dipole distribution, in- 
cluding the 68% spread in the uncertainty of the expecta- 
tion, is shown in Figure 2. This figure shows that the obser- 
vation is consistent with the expectation in both amplitude 
and phase. 
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4 Interference between the Sidereal and So- 
lar Anisotropy 

4.1 The Sidereal Anisotropy Interference in the 
Solar Dipole Anisotropy 

A check applied to ensure that the solar dipole was well 
understood was to compare the expected solar dipole dis- 
tribution to the observed distribution over the three month 
intervals (February-April, May-July, August-October, and 
November- January) with approximately the same detec- 
tor Uvetime. However, complications exist in this simple 
dipole picture because if the data are not collected within 
an integer number of full years, the solar dipole is expected 
to be strongly distorted by the sidereal anisotropy [7]. This 
spurious effect is expected because the sidereal reference 
frame is defined where the celestial sky is fixed, while the 
solar reference frame is defined where the Sun is fixed. It 
takes about four fewer minutes to complete a sidereal day 
than a solar day. Therefore, a static point in the sidereal 
reference frame will move across the solar frame and re- 
turn to the same position on the sky in one full year (this 
is observed implicitly in Figure 2). Therefore, any static 
sidereal distribution averages to zero in the solar reference 
frame after one year but not over partial time intervals of 
the year. 

The points with their statistical errors in Figure 3 show the 
a — OiSun projection of the data collected in the time in- 
terval between February and April. The black Une shows 
the solar dipole expectation from the motion of the Earth 
around the Sun for the same time interval. This illustrates 
the effect of the distortion by the sidereal anisotropy in the 
solar dipole. The plot shows that the observation in the 
solar reference frame is not in agreement with what is ex- 
pected from the solar dipole alone and is strongly distorted 
by the sidereal anisotropy. 

To eliminate such contamination, the experimental side- 
real anisotropy distribution was used to determine how it 
would look like in solar reference frame. Specifically, a nu- 
merical calculation was performed where, every 100/xs, an 
event was generated with a unique UTC time and with right 
ascension from the all-year experimental sidereal distribu- 
tion. The corresponding distribution in a — asun was de- 
termined for the February- April time period and is shown 
in Figure 4. Once the sidereal anisotropy effect in the so- 
lar reference frame is known, it is then subtracted from the 
distribution measured in the solar reference frame. 

Figure 3 shows the distorted cosmic ray distribution in the 
solar reference frame, and Figure 5 shows the distorted dis- 
tribution minus the sidereal interference. The plot shows 
that the data, after correcting for the sidereal effect, is in 
agreement with the solar dipole expectation for the period 
from February to April. This effect was also estimated for 
each of the next three time intervals (May- July, August- 
October, and November- January) using the same method. 
The data was found to be in agreement with the expectation 
within statistical fluctuations. 
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Figure 3: This plot shows the relative intensity of the one- 
dimensional projection of right ascension with respect to 
the Sun. The observed cosmic ray events in the period from 
February to April are shown with statistical uncertainties in 
the black crosses, and the black line is the relative intensity 
of the expected solar dipole for the same time period. 
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Figure 4: This plot shows the relative intensity of the side- 
real anisotropy as it would appear in a — asun projection 
from February to April in absence of the solar dipole. The 
error bars are the statistical uncertainties. 

4.2 The Solar Dipole Anisotropy Interference in 
the Sidereal Anisotropy 

Similarly, the effect of the solar dipole interference in the 
sidereal anisotropy was also checked. This check allows 
for a better estimate of the systematic effect of the sea- 
sonal variation dependence in the sidereal anisotropy of the 
cosmic ray arrival direction [7]. If the relative intensity in 
equatorial coordinates is measured over a full year, then the 
observed sidereal anisotropy is devoid of any distortions by 
the solar dipole as explained previously. However, similar 
to the observation of the solar dipole in a quarter-year, we 
observe (as shown in Figure 6) that the sidereal anisotropy 
measured in the three month periods between February and 
March (in gray points) deviates from the full dataset obser- 
vation (in black Une). 

This deviation is smaller than the solar dipole's distortion 
and is consistent with what we would expect. If the same 
sidereal- to- solar transformation procedure used previously 
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Figure 5: This plot shows the relative intensity of the one- 
dimensional projection in right ascension with respect to 
the Sun. The observed cosmic ray events with their statisti- 
cal uncertainties in the period from February to April minus 
the propagated sidereal anisotropy in the solar reference 
frame for same time period are shown in black crosses. 
The black hne is the relative intensity of the expected solar 
dipole. 



Figure 7: This plot shows the relative intensity of the one- 
dimensional projection of right ascension. The gray crosses 
are the observed cosmic ray events in the sidereal reference 
frame in the period from February to April minus the prop- 
agated solar anisotropy in the sidereal reference frame for 
same time period. The black line is the relative intensity 
of the same sidereal anisotropy except over the full year. 
Both the full year and the quarter-year projections are plot- 
ted with their statistical error bars. 
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Figure 6: This plot shows the relative intensity of the one- 
dimensional projection of right ascension. The gray crosses 
are the observed cosmic ray events in the sidereal reference 
frame in the period from February to April. The black line 
is the relative intensity of the same sidereal anisotropy ex- 
cept over the full year. Both the full year and the quarter- 
year projections are plotted with their statistical error bars. 

is appUed, but for the sidereal anisotropy, then the solar 
dipole contribution from the sidereal observation is elim- 
inated as shown in Figure 7. This effect was also cor- 
rected for the next three time intervals (May- July, August- 
October, and November- January) using the same method. 
The corrected sidereal anisotropy over each of the time 
intervals was also found to be in agreement with that ob- 
served over a full year. 



ration of IceCube from May 2009 to May 2010. The solar 
dipole effect is studied in one full year and in 3 months in- 
tervals. The observed solar dipole in these time intervals 
was found to be consistent with what is expected from the 
motion of the Earth around the Sun in both ampUtude and 
phase. The solar dipole effect is a well understood mea- 
surement that confirms, within statistical uncertainty, the 
reliability of the large-scale sidereal cosmic ray anisotropy 
observations. 
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5 Conclusion 

In this contribution, we presented the observation of the so- 
lar dipole using the data collected by the 59-string configu- 
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Abstract: IceTop, the surface component of the IceCube neutrino observatory, is an air shower array installed at 
the Amundson-Scott South Pole station. It consists of 162 ice Cherenkov detectors, each 90 cm deep with surface 
area 2.7 m^. The tanks detect secondary particles produced by cosmic rays interacting in the atmosphere with 
a counting rate exceeding 1 kHz per detector. The statistical precision enabled by the high counting rate of the 
array enables the study of cosmic ray heliospheric disturbances with unprecedented time resolution in the multi- 
GeV energy regime. Since the detector has each different threshold settings, it is also possible to estimate the energy 
spectrum of particles in such events. We illustrate the performance during a Forbush decrease observed in February, 2011. 

Corresponding authors: Takao Kuwabara^ {takao@bartol.udeLedu), PmlEyensoii^ (evenson@udel.edu) 

^ Bartol Research Institute and Dept. of Physics and Astronomy, University of Delaware, Newark, DE 19716, U.S.A. 

Keywords: IceCube, IceTop, Forbush decrease 



1 Introduction 

IceTop, the surface component of the IceCube neutrino ob- 
servatory, is located at the Amundson-Scott station near the 
geographic South Pole (altitude 2835m). It is an air shower 
array consisting of 162 ice Cherenkov detectors, each 90 
cm deep with surface area 2.7 m^. Signals are recorded by 
photomultiphers in two Digital Optical Modules (DOMs) 
frozen into the ice in each detector. Two DOMs in a detec- 
tor are typically operated at high and low gain settings to 
provide appropriate dynamic range. In this paper we con- 
sider only data from the set of 146 "high gain" DOMs that 
were deployed prior to the 2010-11 austral summer The 
large detector volume coupled with the high altitude and 
nearly zero geomagnetic cutoff at the South Pole yield an 
extremely high counting rate that allows the study of cos- 
mic ray time variations with unprecedented accuracy. 

This analysis uses counting rates from the two discrimina- 
tors in the high gain DOMs. They are termed SPE (Sin- 
gle Photo Electron), and MPE (Multi Photo Electron) for 
historical reasons. The SPE discriminators are actually 
set to thresholds ranging between 1 and 20 photoelectrons 
whereas the MPE threshold are all set near 20 photoelec- 
trons. The counting rates at those SPE and MPE threshold 
settings range from 1 to 10 kHz. IceTop detectors respond 
to several components of the secondaries (muons, neutrons, 
and electomagnetic components [1]) produced by interac- 
tions of >0.6 GeV primary cosmic rays with Earth's at- 
mosphere, and the signal spectrum of the detector contains 
information on the spectrum of the primary particles. By 



observing the intensity variation at different threshold set- 
tings we can estimate the energy spectrum of primary par- 
ticles during the heliospheric events. 

At energies up to ~ 100 GeV, primary galactic cosmic rays 
experience significant variations in response to passing so- 
lar wind disturbances such as interplanetary coronal mass 
ejections (ICMEs). ICMEs consist of large eruptions of 
magnetic field and plasma from Sun. ICMEs accompanied 
by a strong shock often have a depleted region of galactic 
cosmic rays behind the shock and within the CME. When 
Earth enters this depleted region, ground-based cosmic ray 
detectors record a "Forbush Decrease" [2]. During the de- 
crease the cosmic ray intensity fluctuates in response to 
magnetic structures embedded in the ICME. Data with high 
statistical precision from IceTop allow study of these struc- 
tures with unprecedented time resolution, coupled with si- 
multaneous spectral information. In this paper, we illus- 
trate the performance of IceTop during a Forbush decrease 
event observed on 18 February 2011. 

2 Observation of 18 February 2011 Forbush 
Decrease 

The 18 February 2011 Forbush decrease, related to the 
X2.2 solar flare that occurred on 15 February 01:44 UT, 
was the largest Forbush decrease since December 2006. 
A summary of observations of this decrease is shown in 
Figure 1. The four top panels display solar wind speed, 
total component (Bt) and North-South component (Bz) of 
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Figure 1 : Observation of 1 8 February 20 1 1 Forbush Decrease. From top to bottom are, solar wind speed, total component 
(Bt, solid line) and North-South component (Bz, dashed line) of interplanetary magnetic field, proton density, proton 
temperature, and average IceTop MPE scaler rate. 



interplanetary magnetic field (IMF) in GSM coordinates, 
proton density and temperature of solar wind . The data 
plotted are primarily one-minute averaged high resolution 
OMNI data [3] from the WIND spacecraft. When WIND 
data are not available (18 February 20:40 UT to 19 Febru- 
ary 9:50 UT) real-time ACE data [4], shifted by 50 min- 
utes to correct for the time difference between satellite and 
Earth, are substituted. The time of the Storm Sudden Com- 
mencement (SSC) associated with the ICME-driven shock 
is at 18 February 01:30 UT. After the shock, the solar wind 
speed exceeds 700 km/s and IMF magnitude reaches 30 nT. 
Signatures of a magnetic cloud (strong magnetic field with 
smooth rotation, low proton temperature) are seen around 
the period 18 February 05:00 UT - 10:00 UT. 

The lowest panel shows the pressure corrected average 
counting rate of all 146 MPE discriminators. The data in 
the figure is also averaged over one-minute, therefore the 
total count rate from 146 MPE discriminators for each data 
point displayed is ^12 mega count/min. One can see not 
only the rapid decrease and gradual recovery of the For- 



bush decrease but a wealth of fine structure within the de- 
crease. The onset of the Forbush decrease is 18 February 
2:19 UT, which is ^50 min later than the shock passage. 
Such a time lag between the SSC and onset of decrease is 
possibly caused by anisotropy of primary cosmic rays [5]. 
Just before 18 February 5:00 UT, there is a second decrease 
at the same time as the start of the magnetic cloud signa- 
ture. This observation is consistent with the model of a 
classical two-step Forbush decrease [6]. The first step is 
caused by the turbulent sheath region behind the shock and 
the second step is caused by the entry into the enhanced 
and closed magnetic field of magnetic cloud (ejecta). The 
rate reaches minimum when Earth is located at the center 
of the magnetic cloud. 

2.1 Rigidity Dependence of Forbush Decrease 

We can track the evolution of the cosmic ray energy spec- 
trum during the decrease by examining the time behavior as 
a function of discriminator threshold. Figure 2 shows the 
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where Pq = I GV. Then the spectrum near Earth is esti- 
mated as 
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Figure 2: Count rate as a function of discriminator setting 
for the DOMs used in this study. Rates are averaged over 
the interval from 17 February 22:00 UT to 18 February 
02:00, which defines the "base count rate" for the study. 
The line shows the expected values based on a simulation. 



total counting rate as a function of threshold for the time 
interval prior to the event that we have chosen to determine 
the base count rate. 

The energy dependence of the Forbush decrease is shown 
in Figure 3. The size of the decrease for each DOM is cal- 
culated as a percentage deviation from the base count rate. 
Four time periods in the main phase of Forbush decrease 
are plotted in the left panel, while four other time periods 
during the recovery phase are plotted in the right panel. The 
time periods corresponding to each measurement is shown 
in the figure. In the main phase, the first period is at the first 
step of the Forbush decrease, while the remaining three are 
in the second step until the maximum decrease. In the re- 
covery phase, the day by day transition is shown. We con- 
firm that the DOMs with higher count rate (lower threshold 
setting) have a larger decrease, as expected. This rigidity 
dependence is more remarkable in the recovery phase. 

2.2 Estimation of primary spectrum variation by 
the Force-field model 

To quantify the spectral changes during the event we use 
the force-field model [7, 8]. Although the force-field model 
basically describes the long term galactic cosmic ray mod- 
ulation by the solar activity in the heliosphere, it provides 
a simple way to characterize the evolution during this tran- 
sient. 

We define the interstellar spectrum J^^ [l/(m^ sr s MeV)] 
as a function of the particle rigidity P [GV] [9]. 
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P + ^J 1 + iPo/iP + cPW 
where (p [MV] is the so-called modulation parameter. 

The IceTop count rate at different PE thresholds is cal- 
culated by the convolution of the yield function Spe{P) 
[1] with the primary cosmic ray spectrum J{P, (j>) defined 
above, as 



100 



Np,{(b) = J SpeiP)J{P,^)dP 



(3) 



Then, by assuming 00 as the modulation level of the pri- 
mary spectrum before the Forbush decrease the rate de- 
crease during the event is described as a function of (f>. In 
this paper the 0o is roughly set to the 300MV for the normal 
galactic cosmic ray spectrum at this solar minimum [10]. 

The derived decrease size at different time periods is shown 
with curves in Figure 3. Modulation parameters <j> for each 
time period are given in the figure, which reach (f)=4Q6MV 
at the maximum decrease. 

We notice that there are significant differences between the 
force-field model and the observed decrease for the high 
count rate DOMs, especially those set to thresholds less 
than 5 PE, shown as open symbols. Pre-event counting 
rates (Figure 2) for these thresholds are also significantly 
higher than the simulation model that is otherwise a good 
fit. The most likely explanation of this is background due 
to ambient radiation, but we have to this point not made an 
allowance for this in our models. Presence of a constant 
background would clearly reduce the size of the Forbush 
decrease as a percentage, which is how the data in Figure 
3 are displayed. Similar effects were seen in our analysis 
of a solar flare event [11]. When more flare events, with 
varying spectra, are observed we will be able to derive self 
consistent corrections for this background. 

Except for the high count rate DOMs, observed decreases 
in the main phase are well reproduced by the primary spec- 
trum defined by the force-field model. This result suggests 
that it is indeed reasonable to apply the force-field approxi- 
mation from the onset of the transient event [12] through 
the time of maximum decrease. However in the recov- 
ery phase the spectra become progressively less well fit by 
the force-field model. This quantifies the observation that 
higher energy galactic cosmic rays typically recover faster 
than lower energy ones, as observed by [5, 13]. 

3 Summary 

IceTop observed the Forbush decrease event associated 
with the X2.2 solar flare in February, 2011. The high 
counting rate detector makes possible to observe the on- 
set time of the Forbush decrease accurately. The rigidity 
dependence of the Forbush decrease is estimated by the de- 
tector response with the primary spectrum expected from 
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Figure 3: Decrease in counting rate for each DOM at different times during the main phase (left) and recovery phase (right) 
of the Forbush decrease. The decrease, expressed as a percentage of the base count rate, measured from 17 February 22:00 
UT to 18 February 02:00 is plotted as a function of the base count rates. DOMs set to very low threshold (<5 PE) are 
shown as open symbols. Lines are expected decrease size calculated from force-field spectrum at the indicated modulation 
parameter </>. 



the force-field approximation. The observed rigidity de- 
pendence of decrease size is well reproduced by the force- 
field primary spectrum during the main phase. By contrast, 
the strong rigidity dependence in the recovery phase of this 
event is not reproduced by it. 

There are some problems remaining. The expected de- 
crease sizes from the assumed primary spectrum overesti- 
mated the decrease in the high counting rate DOMs. Proba- 
bly the counting rates of these DOMs includes background 
that is not accounted for in the calculation of the yield func- 
tion. Our current yield function is not taking account of the 
effect of the different snow accumulation on the tanks ei- 
ther Further work will be required to refine this. We should 
also note that the current analysis didn't correct for effects 
of the primary cosmic ray anisotropy. 
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